Code Packages:
# Load packages
library(gridExtra) # For using grid.table to save tables as images
## Warning: package 'gridExtra' was built under R version 4.4.3
library(flextable) # Alternative method to save df as image
## Warning: package 'flextable' was built under R version 4.4.3
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ purrr::compose() masks flextable::compose()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(zoo)
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(rlang)
##
## Attaching package: 'rlang'
##
## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
library(nortest) # For ad.test
library(car) # For homogeneity of variance with leveneTest
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(rcompanion) # For scheirerRayHare test
## Warning: package 'rcompanion' was built under R version 4.4.3
library(pastecs) # For stat.desc
## Warning: package 'pastecs' was built under R version 4.4.3
##
## Attaching package: 'pastecs'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## The following object is masked from 'package:tidyr':
##
## extract
library(ggpubr) # For assumption plots
## Warning: package 'ggpubr' was built under R version 4.4.3
##
## Attaching package: 'ggpubr'
##
## The following objects are masked from 'package:flextable':
##
## border, font, rotate
library(conover.test) # For the conover.Iman() test function
library(DescTools) # Contain alternative ConoverTest function
## Warning: package 'DescTools' was built under R version 4.4.3
##
## Attaching package: 'DescTools'
##
## The following object is masked from 'package:car':
##
## Recode
library(rstatix) # tidyverse adapted stat tests and mshapiro_test() functions
## Warning: package 'rstatix' was built under R version 4.4.3
##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filter
library(DHARMa) # For multi-level lm models and shapiro test
## Warning: package 'DHARMa' was built under R version 4.4.3
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(patchwork) # For combining ggplots
## Warning: package 'patchwork' was built under R version 4.4.3
# Read and call data into df
# Primary samples df
ABL90 <- read.csv("Raw_Data/ABL90_Raw.csv")
# River's Parturition metadata
partuition_subcat <- read.csv("Raw_Data/River's_Rockfish_Metadata_Parturition_V7.csv")
# Data sifting: ABL90 dataset
# Step 1: Remove missing info
ABL_set1 <- ABL90 %>%
filter(Patient.ID_edited != "")
# Step 2: Separate Patient.ID into columns of sample types: blood plasma and instant freeze plasma
ABL_set2 <- separate(ABL_set1, 'Patient.ID_edited', into = c("Patient.ID_edited", "Sample_Type"), sep = ",")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 4 rows [489, 619, 1034,
## 1177].
# Step 3:
# Convert Patient.ID_edited column data to character type
ABL_set2$Patient.ID_edited <- as.character(ABL_set2$Patient.ID_edited)
# Step 4 & 5: Sussing out specific sample errors
# Remove insufficient samples
ABL_set3 <- ABL_set2 %>%
filter(!is.na("Type")) %>%
filter(!str_detect(Errors.detected.during.measurement, "Insufficient sample"))
# Remove inhomogeneous samples
ABL_set3 <- ABL_set3 %>%
filter(!is.na("Type")) %>%
filter(!str_detect(Errors.detected.during.measurement, "Inhomogeneous sample"))
# Step 6: Filter for only blood samples
ABL_b_samp <- ABL_set3 %>%
filter(Sample_Type == "b")
Join ABL90 df with parturition_subcat df: Making ABL_merged
# Rename columns to match: Change 'ID' in parturition_subcat to 'Patient.ID_edited' so it is ready to join with ABL90 df
# Please note the 'Treatment' col in parturition_subcat does not match with the finalized 'Ambient.Or.OAH' col in metadata_atresia_guide or ABL90_merged df
# Rename 'ID' col to 'patient.ID_edited'
partuition_subcat <- partuition_subcat %>%
rename(Patient.ID_edited = ID)
# Connect main ABL90 dataset with my (River's) sub categorized parturition metadata
# ABL_merged <- partuition_subcat %>%
# left_join(ABL_b_samp, by = "Patient.ID_edited")
ABL_merged <- ABL_b_samp %>%
inner_join(partuition_subcat, by = "Patient.ID_edited")
Rename Columns Neatly:
# Rename Treatment & Parturition Outcome
ABL_merged <- ABL_merged %>%
rename(Ambient.Or.OAH = Consensus_General_Treatment,
Pregnant.Or.Atresia = Consensus_Atresia_Or_Pregnant)
Remove Replicate ID’s
# Cross Check all Columns: Cross validation to suss correct replicate row
check_params <- ABL_merged %>%
select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
filter(Patient.ID_edited == "9782D")
print(check_params)
## Patient.ID_edited Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1 9782D 12/28/2023 10:04 1016 C 65uL Ambient
## 2 9782D 1/9/2024 17:21 863 S 65uL Ambient
## Pregnant.Or.Atresia pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1 Post Parturition 7.348 0.5 9.8 184 149
## 2 Post Parturition 7.446 0.7 13.3 181 161
## K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1 2.8 1.32 4.2
## 2 2.9 1.31 6.8
check_params <- ABL_merged %>%
select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
filter(Patient.ID_edited == "9783D")
print(check_params)
## Patient.ID_edited Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1 9783D 1/19/2024 10:37 934 S 65uL Ambient
## 2 9783D 1/7/2024 12:32 860 S 65uL Ambient
## Pregnant.Or.Atresia pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1 Post Parturition 7.484 1.7 0.2 175 193
## 2 Post Parturition 7.450 1.3 14.7 179 164
## K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1 34.5 1.29 8.5
## 2 2.5 1.30 5.3
# Removed replicates: 9782D (2x) and 9783D (2x)
ABL_merged <- ABL_merged %>%
filter(Sample.. != "863",
Sample.. != "934")
Review Replicates:
| Row | ID | Time | Sample.. | Measuring.Mode | pH | Glu.mmol.L. | Hct…. | Na…mmol.L. | Cl…mmol.L. | K…mmol.L. | Ca.mmol.L. |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 32 | 9782D | 12/28/2023 10:04 | 1016 | C 65uL | 7.348 | 0.5 | 9.8 | 184 | 149 | 2.8 | 1.32 |
| 39 | 9782D | 1/9/2024 17:21 | 863 | S 65uL | 7.446 | 0.7 | 13.3 | 181 | 161 | 2.9 | 1.31 |
| 36 | 9783D | 1/19/2024 10:37 | 934 | S 65uL | 7.484 | 1.7 | 0.2 | 175 | 193 | 34.5 | 1.29 |
| 41 | 9783D | 1/7/2024 12:32 | 860 | S 65uL | 7.450 | 1.3 | 14.7 | 179 | 164 | 2.5 | 1.30 |
Methods for replicate sample removal: Check processing date of IDs, and remove samples that do not match that date.
Case ID_9782D:
Origianl ABL90 set contains:
9782D,b (Sample.. 1016, Time = 12/28/2023 10:04:00 AM)
9782D,p (Sample.. 1020, Time = 12/28/2023 10:25:00 AM)
9782D,p (Sample.. 865, Time = 1/9/2024 5:29:00 PM)
9782D,b (Sample.. 863, Time = 1/9/2024 5:21:00 PM)
2024_Fish_Metadata.csv says ID_9782D was processed 3/29/2024
Note Processing date in metadata will not match time date stamp in ABL90 data due to the machine being reset, and innacurrate dates being automatically recorded.
Deciding info: - For ID-9782D, refer to metadata sheet and look for other samples with that same processing date, then view thoes samples ABL90 ‘Time’. Try to find a matching ‘Time’ with other IDs and 9782D that were processed together during the same time recorded by the machine.
ID_97838,b (single obs, Sample.. 1024) Time = 12/28/2023 11:49:00 AM
Since 9782D,b Sample.. 1016 and 97838,b Sample.. 1024 both have a date timestamp occurrence on 12/28/2023, I have decided to keep 9782D,b Sample.. 1016 and drop 9783D,b Sample.. 863
Case ID_9783D
Original ABL90 dataset contains:
9783D,b (Sample.. 934, Time = 1/19/2024 10:37:00 AM)
9783D,c (Sample.. 861, Time = 1/7/2024 12:38:00 PM)
9783D,b (Sample.. 860, Time = 1/7/2024 12:32:00 PM)
2024_Fish_Metadata.csv says ID_9783D was processed 3/20/2024
Informed that ID_9783D: Sample.. 934 is a plasma sample (on account of its low hct value), remove this and keep Sample.. 860 (parameter values appear correct).
Remove Motalities and No info IDs and assign analysis ready data-frames:
# New df with Moralities removed: Note none of these samples made it into ABL90 df anyway, so looks like they are already filtered out.
Live_Samples <- ABL_merged %>%
filter(Patient.ID_edited != "9780C", # Mortality
Patient.ID_edited != "777AE", # Mortality
Patient.ID_edited != "777CA") # Mortality after parturition
# New df with mortality and 'No info' Id's removed
Primary_Samples <- Live_Samples %>%
filter(Patient.ID_edited != "777A0", # No info
Patient.ID_edited != "9782F", # No info
Patient.ID_edited != "777B3", # No info
Patient.ID_edited != "777AA", # No info
Patient.ID_edited != "777DE", # No info
Patient.ID_edited != "777CE", # No info
Patient.ID_edited != "777A6") # No info
# New df of Only Ambient Treatment: For testing parturition success
Ambient_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient")
# New df of Fecundity Samples:
Fecundity_Samples <- Primary_Samples %>%
filter(Fecundity_Or_Brood_Count != "NA",
Fecundity_Class != "NA") # removes 97706
# Fecundity samples without atresia Ids
Fecundity_No_Atresia_Samples <- Primary_Samples %>%
filter(All_Fecundity != "Na",
All_Fecundity != 0) # removes atresia IDs
# Fecundity Ambient Samples
Amb_Fec_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient",
All_Fecundity != "NA") # removes OAH treatment samples and ID 97706
# Fecundity Ambient Samples without atresia Ids
Amb_Fec_No_Atresia_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient", # removes OAH treatment
All_Fecundity != "NA", # removes all missing info IDs
All_Fecundity != 0) # removes atresia IDs
Change Data Types:
# Change Continuous Data to type to numeric, double, or factor
# Change data type of ions to factor or numeric
#glimpse(General_Samples)
#glimpse(Ambient_Samples)
Primary_Samples <- Primary_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Ambient_Samples <- Ambient_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Fecundity_Samples <- Fecundity_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Fecundity_No_Atresia_Samples <- Fecundity_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Amb_Fec_Samples <- Amb_Fec_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
# glimpse(Primary_Samples)
# glimpse(Ambient_Samples)
# glimpse(Fecundity_Samples)
# glimpse(Amb_Fec_Samples)
# Fecundity data
#unique(Fecundity_Samples$Fecundity_Class)
# Arrange the order of parturition categories for visualizations & tests
# Primary Samples
# Arrange Treatment
Primary_Samples$Ambient.Or.OAH <- ordered(Primary_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))
# Arrange Parturition outcome
Primary_Samples$Pregnant.Or.Atresia <- ordered(Primary_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))
# Arrange Brood Condition
Primary_Samples$Consensus_Brood_Condition <- ordered(Primary_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Very Low", "Atresia"))
# Arrange Fecundity Class
Primary_Samples$Fecundity_Class <- ordered(Primary_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))
# Remove NAs from ordered character factors
Primary_Samples <- Primary_Samples %>%
filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
#class(Primary_Samples$Fecundity_Class)
# Save merged dataset in a worked folder
write.csv(x = Primary_Samples, file = "Worked-Data/Primary_Samples")
write.csv(x = Ambient_Samples, file = "Worked-Data/Ambient_Samples")
write.csv(x = Fecundity_Samples, file = "Worked-Data/Fecundity_Samples")
write.csv(x = Fecundity_No_Atresia_Samples, file = "Worked-Data/Fecundity_No_Atresia_Samples")
write.csv(x = Amb_Fec_Samples, file = "Worked-Data/Amb-Fec_Samples")
write.csv(x = Amb_Fec_No_Atresia_Samples, file = "Worked-Data/Amb_Fec_No_Atresia_Samples")
# Sample Size
# Primary data: Treatment & Pregnancy outcome
primary_sample_table <- Primary_Samples %>%
group_by(Ambient.Or.OAH, Pregnant.Or.Atresia) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(primary_sample_table)
## # A tibble: 4 × 3
## Ambient.Or.OAH Pregnant.Or.Atresia n_size
## <ord> <ord> <int>
## 1 Ambient Post Parturition 8
## 2 Ambient Atresia 13
## 3 OAH Post Parturition 2
## 4 OAH Atresia 5
pdf("data-images/primary_sample_table.pdf")
grid.table(primary_sample_table)
dev.off()
## png
## 2
png("data-images/primary_sample_table.png")
grid.table(primary_sample_table)
dev.off()
## png
## 2
# Primary data: Brood Condition
primary_sample_table_brood_codnition <- Primary_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(primary_sample_table_brood_codnition)
## # A tibble: 6 × 3
## Ambient.Or.OAH Consensus_Brood_Condition n_size
## <ord> <ord> <int>
## 1 Ambient Excellent 1
## 2 Ambient Normal 4
## 3 Ambient Low 3
## 4 Ambient Atresia 13
## 5 OAH Excellent 2
## 6 OAH Atresia 5
# Primary data: Fecundity Class
primary_sample_table_fecundity_class <- Primary_Samples %>%
group_by(Ambient.Or.OAH, Pregnant.Or.Atresia, Fecundity_Class) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(primary_sample_table_fecundity_class)
## # A tibble: 5 × 4
## Ambient.Or.OAH Pregnant.Or.Atresia Fecundity_Class n_size
## <ord> <ord> <ord> <int>
## 1 Ambient Post Parturition High (>50,000) 4
## 2 Ambient Post Parturition Low (~1,000s) 4
## 3 Ambient Atresia Atresia 13
## 4 OAH Post Parturition High (>50,000) 2
## 5 OAH Atresia Atresia 5
Sample Size Barplot:
# View Sample Distributions
# Sample Size barplot
# Primary Samples: General sample size
primary_sample_barplot <- Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Pregnant.Or.Atresia, y = n_size)) +
geom_bar(stat = "identity" ) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
facet_grid(~ Ambient.Or.OAH) +
labs(title = "Primary Sample Size",
x = "Parturition Status & Treatment Type",
y = "Sample Size (n)") +
guides(fill = guide_legend((title = "Reproductive State"))) +
theme_classic() +
theme(plot.title = element_text(size = 30, hjust = 0.5))
print(primary_sample_barplot)
ggsave(filename = "data-images/.png", plot = primary_sample_barplot, device = "png")
## Saving 7 x 5 in image
ggsave(filename = "data-images/primary_sample_barplot.pdf", plot = primary_sample_barplot, device = "pdf")
## Saving 7 x 5 in image
# Primary samples: Fecundity Class sample size
primary_sample_fecundity_class_barplot <- Primary_Samples %>%
group_by(Fecundity_Class, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Fecundity_Class, y = n_size)) +
geom_bar(stat = "identity" ) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
facet_grid(~ Ambient.Or.OAH) +
labs(title = "Sample Size",
x = "Parturition Status & Treatment Type",
y = "Sample Size (n)") +
guides(fill = guide_legend((title = "Reproductive State"))) +
theme_classic() +
theme(plot.title = element_text(size = 30, hjust = 0.5))
print(primary_sample_fecundity_class_barplot)
# Primary samples: Brood Condition sample size
primary_sample_brood_condition_barplot <- Primary_Samples %>%
group_by(Consensus_Brood_Condition, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Consensus_Brood_Condition, y = n_size)) +
geom_bar(stat = "identity" ) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
facet_grid(~ Ambient.Or.OAH) +
labs(title = "Sample Size",
x = "Parturition Status & Treatment Type",
y = "Sample Size (n)") +
guides(fill = guide_legend((title = "Reproductive State"))) +
theme_classic() +
theme(plot.title = element_text(size = 30, hjust = 0.5))
print(primary_sample_brood_condition_barplot)
pH Summary Stats
# pH
# Summary stats
# Primary data
Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarize(count = n(),
median = round(median(pH), 3),
mean = round(mean(pH), 3),
sd = round(sd(pH), 3),
cv = round(sd(pH)/mean(pH), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.
## # A tibble: 4 × 7
## Pregnant.Or.Atresia Ambient.Or.OAH count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition Ambient 8 7.44 7.44 0.063 0.008
## 2 Post Parturition OAH 2 7.45 7.45 0.057 0.008
## 3 Atresia Ambient 13 7.37 7.38 0.1 0.014
## 4 Atresia OAH 5 7.24 7.25 0.058 0.008
pH Boxplot
# pH boxplot
# Primary samples
pH_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_boxplot)
ggsave(filename = "data-images/pH_primary_boxplot.pdf", plot = pH_primary_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_primary_boxplot.png", plot = pH_primary_boxplot, device = "png")
## Saving 7 x 5 in image
pH Scatterplot: No Atresia IDs
# pH
# Primary Samples
# Weight adjusted scatterplot
# Primary scatterplot: No atresia ids
pH_primary_weight_adj_scatterplot <- ggplot(data = Primary_Samples) +
geom_point(aes(x = pH, y = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight), colour = "black") +
labs(title = "pH", x = "pH", y = "Weight Adjusted Fecundity") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_weight_adj_scatterplot)
## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/pH_primary_weight_adj_scatterplot.pdf", plot = pH_primary_weight_adj_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/pH_primary_weight_adj_scatterplot.png", plot = pH_primary_weight_adj_scatterplot, device = "png")
## Saving 7 x 5 in image
## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Facet by Treatment
pH_primary_weight_adj_scatterplot <- ggplot(data = Primary_Samples) +
geom_point(aes(x = pH, y = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight), colour = "black") +
labs(title = "pH", x = "pH", y = "Weight Adjusted Fecundity") +
scale_y_continuous() +
theme_classic() +
facet_wrap(~ Ambient.Or.OAH) +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_weight_adj_scatterplot)
## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).
Data Distribution
# pH
# Parametric Assumptions
# Data Distribution
hist(Primary_Samples$pH)
plotNormalHistogram(Primary_Samples$pH)
plotNormalDensity(Primary_Samples$pH)
# pH
# Primary Samples
# Interactive aov model
pH_primary_aov_int <- aov(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(pH_primary_aov_int) # interaction almost significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 0.06629 0.06629 9.638 0.00483 **
## Ambient.Or.OAH 1 0.03877 0.03877 5.637 0.02593 *
## Pregnant.Or.Atresia:Ambient.Or.OAH 1 0.02078 0.02078 3.022 0.09495 .
## Residuals 24 0.16506 0.00688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
# pH_primary_aov_additive <- aov(pH ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
# summary(pH_primary_aov_additive)
# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(pH_primary_aov_int)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
##
## $Pregnant.Or.Atresia
## diff lwr upr p adj
## Atresia-Post Parturition -0.1015444 -0.1690501 -0.03403877 0.004834
##
## $Ambient.Or.OAH
## diff lwr upr p adj
## OAH-Ambient -0.08561481 -0.1603143 -0.01091529 0.026428
##
## $`Pregnant.Or.Atresia:Ambient.Or.OAH`
## diff lwr
## Atresia:Ambient-Post Parturition:Ambient -0.06418269 -0.1669825
## Post Parturition:OAH-Post Parturition:Ambient 0.00862500 -0.1722337
## Atresia:OAH-Post Parturition:Ambient -0.19247500 -0.3228940
## Post Parturition:OAH-Atresia:Ambient 0.07280769 -0.1009557
## Atresia:OAH-Atresia:Ambient -0.12829231 -0.2486791
## Atresia:OAH-Post Parturition:OAH -0.20110000 -0.3925028
## upr p adj
## Atresia:Ambient-Post Parturition:Ambient 0.03861711 0.3345551
## Post Parturition:OAH-Post Parturition:Ambient 0.18948366 0.9991631
## Atresia:OAH-Post Parturition:Ambient -0.06205597 0.0023249
## Post Parturition:OAH-Atresia:Ambient 0.24657107 0.6594942
## Atresia:OAH-Atresia:Ambient -0.00790551 0.0337411
## Atresia:OAH-Post Parturition:OAH -0.00969719 0.0369572
# Normality plot
# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
pH_primary_res_qqplot <- ggqqplot(residuals(aov(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Primary pH Interactive Residual QQplot",
subtitle = "residuals(aov(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
x = "pH Theoretical Quantiles (Predicted values)",
y = "pH Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(pH_primary_res_qqplot)
# Normality Test
shapiro.test(Primary_Samples$pH)
##
## Shapiro-Wilk normality test
##
## data: Primary_Samples$pH
## W = 0.97662, p-value = 0.7632
shapiro.test(residuals(lm(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.96397, p-value = 0.431
# Parametric variance test:
bartlett.test(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # nonparametric
##
## Bartlett test of homogeneity of variances
##
## data: pH by Pregnant.Or.Atresia
## Bartlett's K-squared = 3.4356, df = 1, p-value = 0.0638
bartlett.test(pH ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: pH by Ambient.Or.OAH
## Bartlett's K-squared = 0.33993, df = 1, p-value = 0.5599
# Nonparametric variance test:
LeveneTest(pH ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.3324 0.1388
## 26
leveneTest(pH ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1912 0.6655
## 26
# Nonparametric Stat test & Post-Hoc:
kruskal.test(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # Significant
##
## Kruskal-Wallis rank sum test
##
## data: pH by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 6.954, df = 1, p-value = 0.008363
kruskal.test(pH ~ Ambient.Or.OAH, data = Primary_Samples) # Significant
##
## Kruskal-Wallis rank sum test
##
## data: pH by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 3.7502, df = 1, p-value = 0.0528
DunnTest(pH ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Atresia-Post Parturition -8.555556 0.0084 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DunnTest(pH ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## OAH-Ambient -6.952381 0.0528 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Parametric Results: Ambient Fecundity No Atresia Sample
A significant difference detected in pH between pregnancy groups and treatment groups
Interaction almost significant (), interactive model used
Parametric post Hoc: Tukey HSD with interactive model
Ask how to interpret interactive results…
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Met
Homoscedasticity: Not Met, nonparametric
Evidence: data: pH by Pregnant.Or.Atresia Bartlett’s K-squared = 3.4356, df = 1, p-value = 0.0638
Nonparametric variance test (Levene’s test): Met
Non-Parametric Results: - A significant increase detected in pH for post parturition group compared to atresia group. - Evidence: Dunn’s test pval = 0.0084 **
# ph
# Primary Samples
# Regression Model
# Install flexplot
# devtools::install_github("dustinfife/flexplot", ref="development")
require(flexplot)
## Loading required package: flexplot
##
## Attaching package: 'flexplot'
## The following object is masked from 'package:ggplot2':
##
## flip_data
# If interaction not significant, use additive model
# p-value significant if (a < 0.1 or a < 0.05)
pH_primary_lm <- lm(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(pH_primary_lm)
##
## Call:
## lm(formula = pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.149692 -0.042344 -0.008546 0.046671 0.222308
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.37887 0.01969 374.762 < 2e-16 ***
## Pregnant.Or.Atresia.L -0.09379 0.02785 -3.368 0.00255 **
## Ambient.Or.OAH.L -0.04231 0.02785 -1.519 0.14172
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.06846 0.03938 -1.738 0.09495 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08293 on 24 degrees of freedom
## Multiple R-squared: 0.4326, Adjusted R-squared: 0.3617
## F-statistic: 6.099 on 3 and 24 DF, p-value: 0.003094
pH_primary_lm_res_plot <- simulateResiduals(fittedModel = pH_primary_lm)
plot(pH_primary_lm_res_plot)
visualize(pH_primary_lm, plot = "model")
# Residuals of model for shapiro test
pH_primary_lm_res <- residuals(pH_primary_lm)
# Normality test
shapiro.test(pH_primary_lm_res)
##
## Shapiro-Wilk normality test
##
## data: pH_primary_lm_res
## W = 0.96397, p-value = 0.431
# Scedasticity test:
# bartlett's test parametric
bartlett.test(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # nonparametric
##
## Bartlett test of homogeneity of variances
##
## data: pH by Pregnant.Or.Atresia
## Bartlett's K-squared = 3.4356, df = 1, p-value = 0.0638
bartlett.test(pH ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: pH by Ambient.Or.OAH
## Bartlett's K-squared = 0.33993, df = 1, p-value = 0.5599
# Levene's test if nonparametric
leveneTest(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.6518 0.5896
## 24
leveneTest(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # Close to significant
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.3324 0.1388
## 26
leveneTest(pH ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1912 0.6655
## 26
# Non-parametric regression test???:
Results: Primary Sample (With both treatments and parturition groups)
Linear model interaction significant (Pr(>|t|) = 0.09495 .), interactive model used.
Linear Model: Interactive lm param ~ parturition * treatment - A significant relationship detected between pH, parturition groups, and treatment - Evidence: summary lm Call: lm(formula = pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) Residual standard error: 0.08293 on 24 degrees of freedom Multiple R-squared: 0.4326, Adjusted R-squared: 0.3617 F-statistic: 6.099 on 3 and 24 DF, p-value: 0.003094
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality (of residuals) : Met
Homoscedasticity: Met
# ph
# Primary Samples
# Regression Model pt 2
# Individually assess parameter to parturition groups & treatment groups
pH_primary_lm1 <- lm(pH ~ Pregnant.Or.Atresia, data = Primary_Samples)
pH_primary_lm2 <- lm(pH ~ Ambient.Or.OAH, data = Primary_Samples)
# p-value significant if (a < 0.1 or a < 0.05)
summary(pH_primary_lm1) # Significant
##
## Call:
## lm(formula = pH ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.152056 -0.057714 0.000944 0.042900 0.257944
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.39183 0.01833 403.29 <2e-16 ***
## Pregnant.Or.Atresia.L -0.07180 0.02592 -2.77 0.0102 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09295 on 26 degrees of freedom
## Multiple R-squared: 0.2279, Adjusted R-squared: 0.1982
## F-statistic: 7.673 on 1 and 26 DF, p-value: 0.01021
summary(pH_primary_lm2) # Significant
##
## Call:
## lm(formula = pH ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.17414 -0.06893 -0.01614 0.06336 0.19786
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.35350 0.02111 348.403 <2e-16 ***
## Ambient.Or.OAH.L -0.06738 0.02985 -2.257 0.0326 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09672 on 26 degrees of freedom
## Multiple R-squared: 0.1639, Adjusted R-squared: 0.1317
## F-statistic: 5.095 on 1 and 26 DF, p-value: 0.03262
pH_primary_lm_res_plot1 <- simulateResiduals(fittedModel = pH_primary_lm1)
plot(pH_primary_lm_res_plot1)
pH_primary_lm_res_plot2 <- simulateResiduals(fittedModel = pH_primary_lm2)
plot(pH_primary_lm_res_plot2)
# For shapiro test
pH_primary_lm_res1 <- residuals(pH_primary_lm1)
pH_primary_lm_res2 <- residuals(pH_primary_lm2)
# Normality test
shapiro.test(pH_primary_lm_res1) # Normal
##
## Shapiro-Wilk normality test
##
## data: pH_primary_lm_res1
## W = 0.96109, p-value = 0.3699
shapiro.test(pH_primary_lm_res2) # Normal
##
## Shapiro-Wilk normality test
##
## data: pH_primary_lm_res2
## W = 0.97843, p-value = 0.8112
Results: Primary Sample (With both treatments and parturition groups)
Linear Model: Individual test correlation between param and treatment, and param and parturition outcome
Correlation: pH and Post-parturition vs atresia - A significant relationship detected with pH and parturition groups - Evidence: Residual standard error: 0.09295 on 26 degrees of freedom Multiple R-squared: 0.2279, Adjusted R-squared: 0.1982 F-statistic: 7.673 on 1 and 26 DF, p-value: 0.01021
Warning: DHARMA states no significant (n.s.) within group deviation from uniformity, and the Levene Test for homogeneity of variance n.s.
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Met
Homoscedasticity: Met
Residuals Plot (if significant)
# pH
# Primary samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# lm for pH & parturition scatterplot:
pH_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Parturition Group", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_lm_scatterplot1)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/pH_primary_lm_scatterplot1.pdf", plot = pH_primary_lm_scatterplot1, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/pH_primary_lm_scatterplot1.png", plot = pH_primary_lm_scatterplot1, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# lm for pH & parturition scatterplot: faceted by treatment
pH_primary_lm_scatterplot1.1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Parturition Group", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
facet_wrap(~ Ambient.Or.OAH) +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_lm_scatterplot1.1)
## `geom_smooth()` using formula = 'y ~ x'
# lm for pH & treatment scatterplot:
pH_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Treatment Group", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_lm_scatterplot2)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/pH_primary_lm_scatterplot2.pdf", plot = pH_primary_lm_scatterplot2, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/pH_primary_lm_scatterplot2.png", plot = pH_primary_lm_scatterplot2, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# lm for pH & treatment scatterplot: faceted by parturition
pH_primary_lm_scatterplot2.1 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Treatment Group", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
facet_wrap(~ Pregnant.Or.Atresia) +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_lm_scatterplot2.1)
## `geom_smooth()` using formula = 'y ~ x'
# pH
# lm Boxplot
pH_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Parturition Group", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_lm_boxplot1)
## `geom_smooth()` using formula = 'y ~ x'
pH_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = pH)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Treatment Group", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_lm_boxplot2)
## `geom_smooth()` using formula = 'y ~ x'
pH_primary_lm_boxplot_combo <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_primary_lm_boxplot_combo)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
# pH
# lm Boxplot
# Include both regression lines between parturition & treatment groups
ggplot(data = Primary_Samples) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia)) +
geom_point(aes(x = Pregnant.Or.Atresia, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH)) +
geom_point(aes(x = Ambient.Or.OAH, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black")
ggplot(data = Primary_Samples) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia)) +
geom_point(aes(x = Pregnant.Or.Atresia, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH)) +
geom_point(aes(x = Ambient.Or.OAH, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_smooth(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia, group = 1), method = "lm") +
geom_smooth(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH, group = 1), method = "lm")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
ggplot(data = Primary_Samples) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia)) +
geom_point(aes(x = Pregnant.Or.Atresia, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH)) +
geom_point(aes(x = Ambient.Or.OAH, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_smooth(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia, group = 1), method = "lm") +
geom_smooth(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH, group = 1), method = "lm") +
labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
ggplot(data = Primary_Samples) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
geom_point(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Pregnant.Or.Atresia)) +
geom_point(aes(x = Ambient.Or.OAH, y = pH, fill = Pregnant.Or.Atresia), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
geom_smooth(aes(x = Pregnant.Or.Atresia, y = pH, group = 1), method = "lm") +
geom_smooth(aes(x = Ambient.Or.OAH, y = pH, group = 1), method = "lm") +
labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Hct Summary Stats
# hct
# Summary stats
# Primary data
Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarize(count = n(),
median = round(median(Hct....), 3),
mean = round(mean(Hct....), 3),
sd = round(sd(Hct....), 3),
cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.
## # A tibble: 4 × 7
## Pregnant.Or.Atresia Ambient.Or.OAH count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition Ambient 8 14.3 15 3.70 0.247
## 2 Post Parturition OAH 2 16.3 16.3 2.69 0.165
## 3 Atresia Ambient 13 17.3 18.2 2.33 0.128
## 4 Atresia OAH 5 18.6 20 4.50 0.225
Hct Boxplot
# hct boxplot
# Primary samples
hct_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Hct...., fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Hct...., fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
labs(title = "Hematocrit", x = "Parturition Type", y = "Hematocrit (%)") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_primary_boxplot)
ggsave(filename = "data-images/hct_primary_boxplot.pdf", plot = hct_primary_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_primary_boxplot.png", plot = hct_primary_boxplot, device = "png")
## Saving 7 x 5 in image
Data Distribution
# hct
# Parametric Assumptions
# Data Distribution
hist(Primary_Samples$Hct....)
plotNormalHistogram(Primary_Samples$Hct....)
plotNormalDensity(Primary_Samples$Hct....)
Differences:
# hct
# Primary Samples
# Interactive aov model
hct_primary_aov_int <- aov(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_aov_int) # interaction almost significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 76.07 76.07 7.334 0.0123 *
## Ambient.Or.OAH 1 14.13 14.13 1.362 0.2547
## Pregnant.Or.Atresia:Ambient.Or.OAH 1 0.28 0.28 0.027 0.8715
## Residuals 24 248.94 10.37
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
hct_primary_aov_additive <- aov(Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_aov_additive)
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 76.07 76.07 7.631 0.0106 *
## Ambient.Or.OAH 1 14.13 14.13 1.417 0.2451
## Residuals 25 249.22 9.97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(hct_primary_aov_additive)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
##
## $Pregnant.Or.Atresia
## diff lwr upr p adj
## Atresia-Post Parturition 3.44 0.8753289 6.004671 0.0106014
##
## $Ambient.Or.OAH
## diff lwr upr p adj
## OAH-Ambient 1.634286 -1.203694 4.472265 0.2467681
# Normality plot
# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
hct_primary_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Primary Hematocrit Interactive Residual QQplot",
subtitle = "residuals(aov(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
x = "Hct Theoretical Quantiles (Predicted values)",
y = "Hct Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(hct_primary_res_qqplot)
# Normality Test
shapiro.test(Primary_Samples$Hct....)
##
## Shapiro-Wilk normality test
##
## data: Primary_Samples$Hct....
## W = 0.9866, p-value = 0.9687
shapiro.test(residuals(lm(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.96395, p-value = 0.4305
# Parametric variance test:
bartlett.test(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Hct.... by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.16489, df = 1, p-value = 0.6847
bartlett.test(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Hct.... by Ambient.Or.OAH
## Bartlett's K-squared = 0.66339, df = 1, p-value = 0.4154
# Nonparametric variance test:
# LeveneTest(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
# leveneTest(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
# Nonparametric Stat test & Post-Hoc:
# kruskal.test(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
# kruskal.test(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
#
# DunnTest(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
# DunnTest(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
Parametric Results: Ambient Fecundity No Atresia Sample
Interaction not significant (Pr(>F) = 0.8715), additive model used.
Additive summary results:
A significant difference detected in pH between pregnancy groups (Pr(>F) = 0.0106 *), but not treatment groups (Pr(>F) = 0.2451)
Parametric post Hoc: Tukey HSD with additive model
Meaning: Atresia group displayed a significant increase in hematocrit compared to post parturition group.
Evidence: Tukey HSD, Atresia-Post Parturition, p adj = 0.0106014
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Met
Homoscedasticity: Met
Non-Parametric Results: NA
Similarities:
# hct
# Primary Samples
# Regression Model
# If interaction not significant, use additive model
# p-value significant if (a < 0.1 or a < 0.05)
hct_primary_lm_int <- lm(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_lm_int)
##
## Call:
## lm(formula = Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.100 -1.950 -0.900 2.025 6.700
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.3750 0.7647 22.723 <2e-16 ***
## Pregnant.Or.Atresia.L 2.4395 1.0814 2.256 0.0335 *
## Ambient.Or.OAH.L 1.0960 1.0814 1.014 0.3209
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L 0.2500 1.5293 0.163 0.8715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.221 on 24 degrees of freedom
## Multiple R-squared: 0.2666, Adjusted R-squared: 0.1749
## F-statistic: 2.908 on 3 and 24 DF, p-value: 0.05532
hct_primary_lm_add <- lm(Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_lm_add)
##
## Call:
## lm(formula = Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0307 -2.0587 -0.9367 1.8064 6.7693
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.4099 0.7198 24.19 <2e-16 ***
## Pregnant.Or.Atresia.L 2.3419 0.8838 2.65 0.0138 *
## Ambient.Or.OAH.L 1.1642 0.9780 1.19 0.2451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.157 on 25 degrees of freedom
## Multiple R-squared: 0.2657, Adjusted R-squared: 0.207
## F-statistic: 4.524 on 2 and 25 DF, p-value: 0.02104
# hct_primary_lm_res_plot <- simulateResiduals(fittedModel = hct_primary_lm_int)
# plot(hct_primary_lm_res_plot)
hct_primary_lm_res_plot <- simulateResiduals(fittedModel = hct_primary_lm_add)
plot(hct_primary_lm_res_plot)
visualize(hct_primary_lm_add, plot = "model")
# Residuals of model for shapiro test
# hct_primary_lm_int_res <- residuals(hct_primary_lm_int) # does not work
hct_primary_lm_add_res <- residuals(hct_primary_lm_add)
# Normality test
#shapiro.test(hct_primary_lm_int_res) # does not work
shapiro.test(hct_primary_lm_add_res)
##
## Shapiro-Wilk normality test
##
## data: hct_primary_lm_add_res
## W = 0.96235, p-value = 0.3959
# Scedasticity test:
# bartlett's test parametric
bartlett.test(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Hct.... by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.16489, df = 1, p-value = 0.6847
bartlett.test(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Hct.... by Ambient.Or.OAH
## Bartlett's K-squared = 0.66339, df = 1, p-value = 0.4154
# Levene's test if nonparametric
# leveneTest(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
# leveneTest(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
# leveneTest(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
# Non-parametric regression test???:
Results: Primary Sample (With both treatments and parturition groups)
Linear Model: Additive lm param ~ parturition + treatment - A significant relationship detected - Evidence: Multiple R-squared: 0.2657, Adjusted R-squared: 0.207 F-statistic: 4.524 on 2 and 25 DF, p-value: 0.02104
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality (of residuals) : Met
Homoscedasticity: Met
# hct
# Primary Samples
# Regression Model pt 2
# Individually assess parameter to parturition groups & treatment groups
hct_primary_lm1 <- lm(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
hct_primary_lm2 <- lm(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
# p-value significant if (a < 0.1 or a < 0.05)
summary(hct_primary_lm1) # Significant
##
## Call:
## lm(formula = Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.360 -2.025 -0.980 2.075 6.600
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.9800 0.6276 27.055 <2e-16 ***
## Pregnant.Or.Atresia.L 2.4324 0.8876 2.741 0.0109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.183 on 26 degrees of freedom
## Multiple R-squared: 0.2241, Adjusted R-squared: 0.1943
## F-statistic: 7.511 on 1 and 26 DF, p-value: 0.01094
summary(hct_primary_lm2) # not significant
##
## Call:
## lm(formula = Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0810 -2.4810 -0.3119 2.6690 6.3571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.9619 0.7646 23.492 <2e-16 ***
## Ambient.Or.OAH.L 1.3873 1.0813 1.283 0.211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.504 on 26 degrees of freedom
## Multiple R-squared: 0.05954, Adjusted R-squared: 0.02336
## F-statistic: 1.646 on 1 and 26 DF, p-value: 0.2108
hct_primary_lm_res_plot1 <- simulateResiduals(fittedModel = hct_primary_lm1)
plot(hct_primary_lm_res_plot1)
hct_primary_lm_res_plot2 <- simulateResiduals(fittedModel = hct_primary_lm2)
plot(hct_primary_lm_res_plot2)
visualize(hct_primary_lm1, plot = "model")
visualize(hct_primary_lm2, plot = "model")
# For shapiro test
hct_primary_lm_res1 <- residuals(hct_primary_lm1)
hct_primary_lm_res2 <- residuals(hct_primary_lm2)
# Normality test
shapiro.test(hct_primary_lm_res1) # Normal
##
## Shapiro-Wilk normality test
##
## data: hct_primary_lm_res1
## W = 0.95107, p-value = 0.2108
shapiro.test(hct_primary_lm_res2) # Normal
##
## Shapiro-Wilk normality test
##
## data: hct_primary_lm_res2
## W = 0.97352, p-value = 0.6773
Results: Primary Sample Regression Analysis pt2
Linear Model: Individual test correlation between param and treatment, and param and parturition outcome
Correlation: Hct and Post-parturition vs atresia - A significant relationship detected with hematocrit and parturition groups - Evidence: Residual standard error: 3.183 on 26 degrees of freedom Multiple R-squared: 0.2241, Adjusted R-squared: 0.1943 F-statistic: 7.511 on 1 and 26 DF, p-value: 0.01094
Warning: DHARMA states no significant (n.s.) within group deviation from uniformity, and the Levene Test for homogeneity of variance n.s.
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Met
Homoscedasticity: Met (check previous code chunk)
Residuals Plot (if significant)
# pH
# Primary samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# lm for hct & parturition scatterplot:
hct_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Hct....)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_primary_lm_scatterplot1)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/hct_primary_lm_scatterplot1.pdf", plot = hct_primary_lm_scatterplot1, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/hct_primary_lm_scatterplot1.png", plot = hct_primary_lm_scatterplot1, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
hct_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Hct....)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_primary_lm_scatterplot2)
## `geom_smooth()` using formula = 'y ~ x'
# pH
# lm Boxplot
hct_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Hct....)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_primary_lm_boxplot1)
## `geom_smooth()` using formula = 'y ~ x'
hct_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Hct....)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_primary_lm_boxplot2)
## `geom_smooth()` using formula = 'y ~ x'
Glucose Summary Stats
# glucose
# Summary stats
# Primary data
Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarize(count = n(),
median = round(median(Glu..mmol.L.), 3),
mean = round(mean(Glu..mmol.L.), 3),
sd = round(sd(Glu..mmol.L.), 3),
cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.
## # A tibble: 4 × 7
## Pregnant.Or.Atresia Ambient.Or.OAH count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition Ambient 8 1.55 1.61 0.591 0.367
## 2 Post Parturition OAH 2 2.75 2.75 1.63 0.591
## 3 Atresia Ambient 13 2 2.1 0.698 0.332
## 4 Atresia OAH 5 2.1 1.98 0.356 0.18
Glucose Boxplot
# glucose boxplot
# Primary samples
glucose_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L., fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L., fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
labs(title = "Glucose", x = "Parturition Type", y = "Glucose (mmol/L)") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_primary_boxplot)
ggsave(filename = "data-images/glucose_primary_boxplot.pdf", plot = glucose_primary_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_primary_boxplot.png", plot = glucose_primary_boxplot, device = "png")
## Saving 7 x 5 in image
Data Distribution
# glucose
# Parametric Assumptions
# Data Distribution
hist(Primary_Samples$Glu..mmol.L.)
plotNormalHistogram(Primary_Samples$Glu..mmol.L.)
plotNormalDensity(Primary_Samples$Glu..mmol.L.)
Differences:
# glucose
# Primary Samples
# Interactive aov model
glucose_primary_aov_int <- aov(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(glucose_primary_aov_int) # interaction significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 0.330 0.3303 0.693 0.4134
## Ambient.Or.OAH 1 0.369 0.3690 0.774 0.3877
## Pregnant.Or.Atresia:Ambient.Or.OAH 1 1.753 1.7533 3.678 0.0671 .
## Residuals 24 11.442 0.4767
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
# glucose_primary_aov_additive <- aov(Glu..mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
# summary(glucose_primary_aov_additive)
# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(glucose_primary_aov_int)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
##
## $Pregnant.Or.Atresia
## diff lwr upr p adj
## Atresia-Post Parturition 0.2266667 -0.3353791 0.7887124 0.4134193
##
## $Ambient.Or.OAH
## diff lwr upr p adj
## OAH-Ambient 0.264127 -0.3578141 0.8860681 0.3894475
##
## $`Pregnant.Or.Atresia:Ambient.Or.OAH`
## diff lwr upr
## Atresia:Ambient-Post Parturition:Ambient 0.4875 -0.3684013 1.3434013
## Post Parturition:OAH-Post Parturition:Ambient 1.1375 -0.3683119 2.6433119
## Atresia:OAH-Post Parturition:Ambient 0.3675 -0.7183564 1.4533564
## Post Parturition:OAH-Atresia:Ambient 0.6500 -0.7967373 2.0967373
## Atresia:OAH-Atresia:Ambient -0.1200 -1.1223290 0.8823290
## Atresia:OAH-Post Parturition:OAH -0.7700 -2.3636015 0.8236015
## p adj
## Atresia:Ambient-Post Parturition:Ambient 0.4131342
## Post Parturition:OAH-Post Parturition:Ambient 0.1869587
## Atresia:OAH-Post Parturition:Ambient 0.7872361
## Post Parturition:OAH-Atresia:Ambient 0.6087444
## Atresia:OAH-Atresia:Ambient 0.9872576
## Atresia:OAH-Post Parturition:OAH 0.5518182
# Normality plot
# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
glucose_primary_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Glucose Interactive Residual QQplot",
subtitle = "residuals(aov(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
x = "Glucose Theoretical Quantiles (Predicted values)",
y = "Glucose Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(glucose_primary_res_qqplot)
# Normality Test
shapiro.test(Primary_Samples$Glu..mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Primary_Samples$Glu..mmol.L.
## W = 0.84841, p-value = 0.0008683
shapiro.test(residuals(lm(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.88835, p-value = 0.00617
# Parametric variance test:
bartlett.test(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Glu..mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 1.6702, df = 1, p-value = 0.1962
bartlett.test(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Glu..mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 0.27111, df = 1, p-value = 0.6026
# Nonparametric variance test:
LeveneTest(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.9338 0.3428
## 26
leveneTest(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0651 0.8005
## 26
# Nonparametric Stat test & Post-Hoc:
kruskal.test(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Glu..mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 2.9147, df = 1, p-value = 0.08778
kruskal.test(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Glu..mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 0.91756, df = 1, p-value = 0.3381
DunnTest(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Atresia-Post Parturition 5.522222 0.0878 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DunnTest(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## OAH-Ambient 3.428571 0.3381
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Parametric Results: Ambient Fecundity No Atresia Sample
Interaction significant (Pr(>F) = 0.0671), interactive model used.
Interactive summary results:
ask?
Parametric post Hoc: Tukey HSD with interactive model
Meaning: No groups appear significantly different from eachother…
Evidence:
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Not met, W = 0.84841, p-value = 0.0008683
Homoscedasticity: Met
Non-Parametric Results:
Kruskal-Wallis rank sum test: - A significant difference in glucose detected between pregnant and atresia groups - Evidence: Kruskal-Wallis chi-squared = 2.9147, df = 1, p-value = 0.08778
Nonparametric Post-Hoc: Dunn’s test - A significant increase in glucose detected in Atresia group compared to post parturition group (mean.rank.diff = 5.522222, pval = 0.0878).
Similarities:
# glucose
# Primary Samples
# Regression Model
# If interaction not significant, use additive model
# p-value significant if (a < 0.1 or a < 0.05)
glucose_primary_lm_int <- lm(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(glucose_primary_lm_int) # interaction significant
##
## Call:
## lm(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1500 -0.3344 -0.0400 0.2050 2.1000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.11062 0.16393 12.875 2.87e-12 ***
## Pregnant.Or.Atresia.L -0.09988 0.23184 -0.431 0.6704
## Ambient.Or.OAH.L 0.35974 0.23184 1.552 0.1338
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.62875 0.32787 -1.918 0.0671 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6905 on 24 degrees of freedom
## Multiple R-squared: 0.1765, Adjusted R-squared: 0.07358
## F-statistic: 1.715 on 3 and 24 DF, p-value: 0.1907
glucose_primary_lm_add <- lm(Glu..mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(glucose_primary_lm_add)
##
## Call:
## lm(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8868 -0.4033 -0.1228 0.1087 2.2073
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0228 0.1656 12.213 4.92e-12 ***
## Pregnant.Or.Atresia.L 0.1456 0.2034 0.716 0.481
## Ambient.Or.OAH.L 0.1882 0.2250 0.836 0.411
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7265 on 25 degrees of freedom
## Multiple R-squared: 0.05033, Adjusted R-squared: -0.02565
## F-statistic: 0.6624 on 2 and 25 DF, p-value: 0.5244
# Simulate residuals plot
# glucose_primary_lm_res_plot <- simulateResiduals(fittedModel = glucose_primary_lm_int)
# plot(glucose_primary_lm_res_plot)
glucose_primary_lm_res_plot <- simulateResiduals(fittedModel = glucose_primary_lm_int)
plot(glucose_primary_lm_res_plot)
visualize(glucose_primary_lm_int, plot = "model")
# Residuals of model for shapiro test
glucose_primary_lm_int_res <- residuals(glucose_primary_lm_int)
glucose_primary_lm_add_res <- residuals(glucose_primary_lm_add)
# Normality test
shapiro.test(glucose_primary_lm_int_res) # works
##
## Shapiro-Wilk normality test
##
## data: glucose_primary_lm_int_res
## W = 0.88835, p-value = 0.00617
shapiro.test(glucose_primary_lm_add_res) # works
##
## Shapiro-Wilk normality test
##
## data: glucose_primary_lm_add_res
## W = 0.81579, p-value = 0.0002053
# Scedasticity test:
# bartlett's test parametric
bartlett.test(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Glu..mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 1.6702, df = 1, p-value = 0.1962
bartlett.test(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Glu..mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 0.27111, df = 1, p-value = 0.6026
# Levene's test if nonparametric
leveneTest(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.9355 0.1508
## 24
leveneTest(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.9338 0.3428
## 26
leveneTest(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0651 0.8005
## 26
# Non-parametric regression test???:
Results: Primary Sample (With both treatments and parturition groups)
Linear Model: Interactive lm param ~ parturition * treatment - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.1765, Adjusted R-squared: 0.07358 F-statistic: 1.715 on 3 and 24 DF, p-value: 0.1907
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality (of residuals) : Not met
Homoscedasticity: Met
Non-parametric Tests:
nonparametric normality test (Levene’s test): Met
Nonparametric regression stat test???
# glucose
# Primary Samples
# Regression Model pt 2
# Individually assess parameter to parturition groups & treatment groups
glucose_primary_lm1 <- lm(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
glucose_primary_lm2 <- lm(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
# p-value significant if (a < 0.1 or a < 0.05)
summary(glucose_primary_lm1) # Significant
##
## Call:
## lm(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9400 -0.3850 -0.1033 0.1333 2.1333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9533 0.1424 13.714 2.06e-13 ***
## Pregnant.Or.Atresia.L 0.1603 0.2014 0.796 0.433
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7223 on 26 degrees of freedom
## Multiple R-squared: 0.02377, Adjusted R-squared: -0.01378
## F-statistic: 0.6331 on 1 and 26 DF, p-value: 0.4334
summary(glucose_primary_lm2) # not significant
##
## Call:
## lm(formula = Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0143 -0.4393 -0.1143 0.1857 2.2857
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0571 0.1570 13.10 5.87e-13 ***
## Ambient.Or.OAH.L 0.2020 0.2221 0.91 0.371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 26 degrees of freedom
## Multiple R-squared: 0.03085, Adjusted R-squared: -0.00643
## F-statistic: 0.8275 on 1 and 26 DF, p-value: 0.3714
glucose_primary_lm_res_plot1 <- simulateResiduals(fittedModel = glucose_primary_lm1)
plot(glucose_primary_lm_res_plot1) # Warning appears
glucose_primary_lm_res_plot2 <- simulateResiduals(fittedModel = glucose_primary_lm2)
plot(glucose_primary_lm_res_plot2)
visualize(glucose_primary_lm1, plot = "model")
visualize(glucose_primary_lm2, plot = "model")
# For shapiro test
glucose_primary_lm_res1 <- residuals(glucose_primary_lm1)
glucose_primary_lm_res2 <- residuals(glucose_primary_lm2)
# Normality test
shapiro.test(glucose_primary_lm_res1) # Non-parametric
##
## Shapiro-Wilk normality test
##
## data: glucose_primary_lm_res1
## W = 0.79757, p-value = 9.675e-05
shapiro.test(glucose_primary_lm_res2) # Non-parametric
##
## Shapiro-Wilk normality test
##
## data: glucose_primary_lm_res2
## W = 0.8557, p-value = 0.001221
Results: Primary Sample Regression Analysis pt2
Linear Model: Individual test correlation between param and treatment, and param and parturition outcome
Correlation: - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: - Meaning:
Warning: DHARMA states for Glu..mmol.L. ~ Pregnant.Or.Atresia QQ plot residuals, KS Test: P = 0.0271 Deviation Significant
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Non-parametric See previous code block
Homoscedasticity: Met see previous code block
Residuals Plot (if significant)
# glucose
# Primary samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# lm for glucose & parturition scatterplot:
glucose_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Parturition Group", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_primary_lm_scatterplot1)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/glucose_primary_lm_scatterplot1.pdf", plot = glucose_primary_lm_scatterplot1, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/glucose_primary_lm_scatterplot1.png", plot = glucose_primary_lm_scatterplot1, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# lm for glucose & treatment scatterplot:
glucose_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Glu..mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Treatment Group", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_primary_lm_scatterplot2)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/glucose_primary_lm_scatterplot2.pdf", plot = glucose_primary_lm_scatterplot2, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/glucose_primary_lm_scatterplot2.png", plot = glucose_primary_lm_scatterplot2, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# glucose
# lm Boxplot
glucose_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Parturition Group", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_primary_lm_boxplot1)
## `geom_smooth()` using formula = 'y ~ x'
glucose_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Glu..mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Parturition Group", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_primary_lm_boxplot2)
## `geom_smooth()` using formula = 'y ~ x'
Sodium Summary Stats
# Na+
# Summary stats
# Primary data
Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarize(count = n(),
median = round(median(Na...mmol.L.), 3),
mean = round(mean(Na...mmol.L.), 3),
sd = round(sd(Na...mmol.L.), 3),
cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.
## # A tibble: 4 × 7
## Pregnant.Or.Atresia Ambient.Or.OAH count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition Ambient 8 176. 176. 3.41 0.019
## 2 Post Parturition OAH 2 178 178 0 0
## 3 Atresia Ambient 13 178 181. 11.4 0.063
## 4 Atresia OAH 5 174 174 2.24 0.013
Sodium Boxplot
# Na+ boxplot
# Primary samples
sodium_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Na...mmol.L., fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Na...mmol.L., fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
labs(title = "Sodium", x = "Parturition Type", y = "Sodium (mmol/L)") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_primary_boxplot)
ggsave(filename = "data-images/sodium_primary_boxplot.pdf", plot = sodium_primary_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_primary_boxplot.png", plot = sodium_primary_boxplot, device = "png")
## Saving 7 x 5 in image
Data Distribution
# Na+
# Parametric Assumptions
# Data Distribution
hist(Primary_Samples$Na...mmol.L.)
plotNormalHistogram(Primary_Samples$Na...mmol.L.)
plotNormalDensity(Primary_Samples$Na...mmol.L.)
Differences:
# Na+
# Primary Samples
# Interactive aov model
sodium_primary_aov_int <- aov(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(sodium_primary_aov_int) # interaction not significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 58.7 58.72 0.849 0.366
## Ambient.Or.OAH 1 97.2 97.24 1.406 0.247
## Pregnant.Or.Atresia:Ambient.Or.OAH 1 99.7 99.66 1.441 0.242
## Residuals 24 1659.8 69.16
# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
sodium_primary_aov_add <- aov(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(sodium_primary_aov_add)
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 58.7 58.72 0.834 0.370
## Ambient.Or.OAH 1 97.2 97.24 1.382 0.251
## Residuals 25 1759.5 70.38
# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(sodium_primary_aov_add)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
##
## $Pregnant.Or.Atresia
## diff lwr upr p adj
## Atresia-Post Parturition 3.022222 -3.792266 9.83671 0.3697553
##
## $Ambient.Or.OAH
## diff lwr upr p adj
## OAH-Ambient -4.287831 -11.82852 3.252855 0.2525973
# Normality plot
# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
sodium_primary_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Sodium Interactive Residual QQplot",
subtitle = "residuals(aov(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
x = "Sodium Theoretical Quantiles (Predicted values)",
y = "Sodium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(sodium_primary_res_qqplot)
# Normality Test
shapiro.test(Primary_Samples$Na...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Primary_Samples$Na...mmol.L.
## W = 0.52152, p-value = 1.916e-08
shapiro.test(residuals(lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.56167, p-value = 5.246e-08
# Parametric variance test:
bartlett.test(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Na...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 10.874, df = 1, p-value = 0.0009751
bartlett.test(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Na...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 8.4559, df = 1, p-value = 0.003639
# Nonparametric variance test:
LeveneTest(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.4826 0.4934
## 26
leveneTest(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2851 0.5979
## 26
# Nonparametric Stat test & Post-Hoc:
kruskal.test(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Na...mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 0.084351, df = 1, p-value = 0.7715
kruskal.test(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Na...mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 2.8469, df = 1, p-value = 0.09155
DunnTest(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Atresia-Post Parturition 0.9333333 0.7715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DunnTest(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## OAH-Ambient -6 0.0916 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sodium_primary_aov_res_plot <- simulateResiduals(fittedModel = sodium_primary_aov_add)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(sodium_primary_aov_res_plot) # KS significant
visualize(sodium_primary_aov_add)
boxplot(data = Primary_Samples, Na...mmol.L. ~ Ambient.Or.OAH)
# Test outliers
outlierTest(sodium_primary_aov_add)
## rstudent unadjusted p-value Bonferroni p
## 24 12.31196 7.3338e-12 2.0535e-10
# Used Ranked model to adjust for outlier
boxplot(rank(Na...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)
boxplot(rank(Na...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)
# Nonparametric Wilcoxon test
wilcox.test(rank(Na...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(Na...mmol.L.) by Pregnant.Or.Atresia
## W = 84, p-value = 0.7901
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(rank(Na...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(Na...mmol.L.) by Ambient.Or.OAH
## W = 105, p-value = 0.09682
## alternative hypothesis: true location shift is not equal to 0
Parametric Results: Ambient Fecundity No Atresia Sample
Warning: Severe outlier may be interfering with results
Interaction not significant (Pr(>F) = 0.242), additive model used.
Additive summary results:
No significant differences between parturition and treatment groups.
Parametric post Hoc: Tukey HSD with interactive model
Meaning: No groups appear significantly different from eachother…
Evidence:
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Really not met, W = 0.52152, p-value = 1.916e-08
Homoscedasticity: Not met (for both factor groups)
Non-Parametric Results:
Kruskal-Wallis rank sum test: - A significant difference in sodium detected between Treatments: Ambient and OAH groups - Evidence: Kruskal-Wallis chi-squared = 2.8469, df = 1, p-value = 0.09155
Nonparametric Post-Hoc: Dunn’s test - A significant difference in sodium detected in treatment groups (mean.rank.diff = -6 pval = 0.0916). - How to know which group is greater? - Generate boxplot for sodium just with treatment groups?
Similarities:
# sodium
# Primary Samples
# Regression Model
# p-value significant if (a < 0.1 or a < 0.05)
sodium_primary_lm_int <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(sodium_primary_lm_int)
##
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.231 -3.481 -1.115 1.063 36.769
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 177.2452 1.9745 89.769 <2e-16 ***
## Pregnant.Or.Atresia.L 0.5235 2.7923 0.187 0.853
## Ambient.Or.OAH.L -1.7610 2.7923 -0.631 0.534
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -4.7404 3.9489 -1.200 0.242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.316 on 24 degrees of freedom
## Multiple R-squared: 0.1335, Adjusted R-squared: 0.02514
## F-statistic: 1.232 on 3 and 24 DF, p-value: 0.3198
sodium_primary_lm_add <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(sodium_primary_lm_add)
##
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.064 -3.422 -1.743 0.907 37.578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 176.583 1.913 92.330 <2e-16 ***
## Pregnant.Or.Atresia.L 2.375 2.348 1.011 0.322
## Ambient.Or.OAH.L -3.055 2.599 -1.175 0.251
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.389 on 25 degrees of freedom
## Multiple R-squared: 0.08142, Adjusted R-squared: 0.007938
## F-statistic: 1.108 on 2 and 25 DF, p-value: 0.3459
# sodium_primary_lm_res_plot <- simulateResiduals(fittedModel = sodium_primary_lm_int)
# plot(sodium_primary_lm_res_plot)
sodium_primary_lm_res_plot <- simulateResiduals(fittedModel = sodium_primary_lm_add)
plot(sodium_primary_lm_res_plot) # KS significant
visualize(sodium_primary_lm_add, plot = "model")
# Check Dispersion or Outliers
plotResiduals(sodium_primary_lm_add)
testDispersion(sodium_primary_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# Transform data: square root
sodium_primary_sqrt <- Primary_Samples %>%
mutate(Na...mmol.L. = sqrt(Na...mmol.L.))
# Rerun lm model with square root transformed data
sodium_primary_sqrt_lm_int <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = sodium_primary_sqrt)
summary(sodium_primary_sqrt_lm_int) # still not significant
##
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = sodium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26562 -0.12425 -0.03893 0.04045 1.30830
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.31135 0.07073 188.202 <2e-16 ***
## Pregnant.Or.Atresia.L 0.01733 0.10003 0.173 0.864
## Ambient.Or.OAH.L -0.06389 0.10003 -0.639 0.529
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.17549 0.14146 -1.241 0.227
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2979 on 24 degrees of freedom
## Multiple R-squared: 0.1388, Adjusted R-squared: 0.03118
## F-statistic: 1.29 on 3 and 24 DF, p-value: 0.3006
sodium_primary_sqrt_lm_add <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = sodium_primary_sqrt)
summary(sodium_primary_sqrt_lm_add)
##
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = sodium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30517 -0.12246 -0.06196 0.03584 1.33823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.28684 0.06864 193.567 <2e-16 ***
## Pregnant.Or.Atresia.L 0.08586 0.08428 1.019 0.318
## Ambient.Or.OAH.L -0.11178 0.09327 -1.198 0.242
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3011 on 25 degrees of freedom
## Multiple R-squared: 0.08361, Adjusted R-squared: 0.01029
## F-statistic: 1.14 on 2 and 25 DF, p-value: 0.3358
# Check hows it
sodium_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = sodium_primary_sqrt_lm_add)
plot(sodium_primary_sqrt_lm_res_plot) # KS still significant
## qu = 0.5, log(sigma) = -3.042338 : outer Newton did not converge fully.
plotResiduals(sodium_primary_sqrt_lm_add)
## qu = 0.5, log(sigma) = -3.042338 : outer Newton did not converge fully.
testDispersion(sodium_primary_sqrt_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
sodium_primary_lm_int_res <- residuals(sodium_primary_lm_int)
sodium_primary_lm_add_res <- residuals(sodium_primary_lm_add)
# Normality test
shapiro.test(sodium_primary_lm_int_res) # significant
##
## Shapiro-Wilk normality test
##
## data: sodium_primary_lm_int_res
## W = 0.56167, p-value = 5.246e-08
shapiro.test(sodium_primary_lm_add_res) # significant
##
## Shapiro-Wilk normality test
##
## data: sodium_primary_lm_add_res
## W = 0.56739, p-value = 6.085e-08
# Scedasticity test:
# bartlett's test parametric
bartlett.test(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Na...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 10.874, df = 1, p-value = 0.0009751
bartlett.test(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Na...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 8.4559, df = 1, p-value = 0.003639
# Levene's test if nonparametric
leveneTest(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.3841 0.7654
## 24
leveneTest(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.4826 0.4934
## 26
leveneTest(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2851 0.5979
## 26
# Non-parametric regression test???:
Results: Primary Sample (With both treatments and parturition groups)
Linear Model: Additive lm param ~ parturition + treatment (no sqrt adjust) - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.08142, Adjusted R-squared: 0.007938 F-statistic: 1.108 on 2 and 25 DF, p-value: 0.3459
Warning: DHARMA states for Na…mmol.L ~ Pregnant.Or.Atresia QQ plot residuals, KS Test: P = 0.0271 Deviation Significant
Correlation: Using sqrt transformed sodium values (additive model) - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: Multiple R-squared: 0.08361, Adjusted R-squared: 0.01029 F-statistic: 1.14 on 2 and 25 DF, p-value: 0.3358
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality (of residuals) : Not met
Homoscedasticity: Not Met
Non-parametric Tests:
nonparametric normality test (Levene’s test): Met
Nonparametric regression stat test???
# glucose
# Primary Samples
# Regression Model pt 2
# Individually assess parameter to parturition groups & treatment groups
sodium_primary_lm1 <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia, data = sodium_primary_sqrt)
sodium_primary_lm2 <- lm(Na...mmol.L. ~ Ambient.Or.OAH, data = sodium_primary_sqrt)
# p-value significant if (a < 0.1 or a < 0.05)
summary(sodium_primary_lm1)
##
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia, data = sodium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30599 -0.11729 -0.04291 0.04232 1.38214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.32812 0.05987 222.606 <2e-16 ***
## Pregnant.Or.Atresia.L 0.07717 0.08467 0.911 0.37
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3036 on 26 degrees of freedom
## Multiple R-squared: 0.03096, Adjusted R-squared: -0.006314
## F-statistic: 0.8306 on 1 and 26 DF, p-value: 0.3705
summary(sodium_primary_lm2)
##
## Call:
## lm(formula = Na...mmol.L. ~ Ambient.Or.OAH, data = sodium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38034 -0.08912 -0.04079 -0.00125 1.38449
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.30708 0.06575 202.382 <2e-16 ***
## Ambient.Or.OAH.L -0.10360 0.09299 -1.114 0.275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3013 on 26 degrees of freedom
## Multiple R-squared: 0.04557, Adjusted R-squared: 0.008856
## F-statistic: 1.241 on 1 and 26 DF, p-value: 0.2754
sodium_primary_lm_res_plot1 <- simulateResiduals(fittedModel = sodium_primary_lm1)
plot(sodium_primary_lm_res_plot1) # Warning appears
sodium_primary_lm_res_plot2 <- simulateResiduals(fittedModel = sodium_primary_lm2)
plot(sodium_primary_lm_res_plot2)
visualize(sodium_primary_lm1, plot = "model")
visualize(sodium_primary_lm2, plot = "model")
# For shapiro test
sodium_primary_lm_res1 <- residuals(sodium_primary_lm1)
sodium_primary_lm_res2 <- residuals(sodium_primary_lm2)
# Normality test
shapiro.test(sodium_primary_lm_res1) # Non-parametric
##
## Shapiro-Wilk normality test
##
## data: sodium_primary_lm_res1
## W = 0.58314, p-value = 9.212e-08
shapiro.test(sodium_primary_lm_res2) # Non-parametric
##
## Shapiro-Wilk normality test
##
## data: sodium_primary_lm_res2
## W = 0.55298, p-value = 4.198e-08
Results: Primary Sample Regression Analysis pt2
Linear Model: Individual test correlation between param and treatment, and param and parturition outcome
Correlation: Using sqrt transformed sodium values - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: - Meaning:
Warning: DHARMA states for Na…mmol.L. ~ Pregnant.Or.Atresia QQ plot residuals, KS Test: P < 0.05. Deviation Significant
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Non-parametric See previous code block
Homoscedasticity: Not Met see previous code block
Residuals Plot (if significant)
# sodium
# Primary samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# lm for sodium & parturition scatterplot:
sodium_primary_lm_scatterplot1 <- ggplot(data = sodium_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Na...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Parturition Group", y = "Sqrt Sodium", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_primary_lm_scatterplot1)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/sodium_primary_lm_scatterplot1.pdf", plot = sodium_primary_lm_scatterplot1, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/sodium_primary_lm_scatterplot1.png", plot = sodium_primary_lm_scatterplot1, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# lm for sodium & treatment scatterplot:
sodium_primary_lm_scatterplot2 <- ggplot(data = sodium_primary_sqrt, aes(x = Ambient.Or.OAH, y = Na...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Treatment Group", y = "Sqrt Sodium", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_primary_lm_scatterplot2)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/sodium_primary_lm_scatterplot2.pdf", plot = sodium_primary_lm_scatterplot2, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/sodium_primary_lm_scatterplot2.png", plot = sodium_primary_lm_scatterplot2, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# sodium
# lm Boxplot
sodium_primary_lm_boxplot1 <- ggplot(data = sodium_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Na...mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Parturition Group", y = "Sqrt Sodium", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_primary_lm_boxplot1)
## `geom_smooth()` using formula = 'y ~ x'
sodium_primary_lm_boxplot2 <- ggplot(data = sodium_primary_sqrt, aes(x = Ambient.Or.OAH, y = Na...mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Parturition Group", y = "Sqrt Sodium", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_primary_lm_boxplot2)
## `geom_smooth()` using formula = 'y ~ x'
Chloride Summary Stats
# Cl-
# Summary stats
# Primary data
Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarize(count = n(),
median = round(median(Cl...mmol.L.), 3),
mean = round(mean(Cl...mmol.L.), 3),
sd = round(sd(Cl...mmol.L.), 3),
cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.
## # A tibble: 4 × 7
## Pregnant.Or.Atresia Ambient.Or.OAH count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition Ambient 8 157 158. 3.10 0.02
## 2 Post Parturition OAH 2 160 160 2.83 0.018
## 3 Atresia Ambient 13 158 158. 4.61 0.029
## 4 Atresia OAH 5 159 159. 1.79 0.011
Chloride Boxplot
# Cl- boxplot
# Primary samples
chloride_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L., fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L., fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
labs(title = "Chloride", x = "Parturition Type", y = "Chloride (mmol/L)") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_primary_boxplot)
ggsave(filename = "data-images/chloride_primary_boxplot.pdf", plot = chloride_primary_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_primary_boxplot.png", plot = chloride_primary_boxplot, device = "png")
## Saving 7 x 5 in image
Data Distribution
# Cl-
# Parametric Assumptions
# Data Distribution
hist(Primary_Samples$Cl...mmol.L.)
plotNormalHistogram(Primary_Samples$Cl...mmol.L.)
plotNormalDensity(Primary_Samples$Cl...mmol.L.)
Differences:
# Cl-
# Primary Samples
# Interactive aov model
chloride_primary_aov_int <- aov(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(chloride_primary_aov_int) # interaction not significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 1.4 1.400 0.098 0.757
## Ambient.Or.OAH 1 7.5 7.536 0.526 0.475
## Pregnant.Or.Atresia:Ambient.Or.OAH 1 2.5 2.533 0.177 0.678
## Residuals 24 343.5 14.314
# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
chloride_primary_aov_add <- aov(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(chloride_primary_aov_add)
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 1.4 1.400 0.101 0.753
## Ambient.Or.OAH 1 7.5 7.536 0.544 0.467
## Residuals 25 346.1 13.843
# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(chloride_primary_aov_add)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
##
## $Pregnant.Or.Atresia
## diff lwr upr p adj
## Atresia-Post Parturition 0.4666667 -2.555517 3.48885 0.7531103
##
## $Ambient.Or.OAH
## diff lwr upr p adj
## OAH-Ambient 1.193651 -2.150597 4.537899 0.4691166
# Normality plot
# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
chloride_primary_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Chloride Interactive Residual QQplot",
subtitle = "residuals(aov(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
x = "Chloride Theoretical Quantiles (Predicted values)",
y = "Chloride Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5))
print(chloride_primary_res_qqplot)
# Normality Test
shapiro.test(Primary_Samples$Cl...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Primary_Samples$Cl...mmol.L.
## W = 0.97188, p-value = 0.6318
shapiro.test(residuals(lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.96898, p-value = 0.5534
# Parametric variance test:
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.76057, df = 1, p-value = 0.3832
bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6573, df = 1, p-value = 0.05582
# Nonparametric variance test:
LeveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1495 0.7021
## 26
leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.9761 0.09637 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Nonparametric Stat test & Post-Hoc:
kruskal.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Cl...mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 0.23333, df = 1, p-value = 0.6291
kruskal.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Cl...mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 1.1429, df = 1, p-value = 0.285
DunnTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Atresia-Post Parturition 1.555556 0.6291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DunnTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## OAH-Ambient 3.809524 0.2850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
visualize(chloride_primary_aov_add)
# Test outliers
outlierTest(chloride_primary_aov_add)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 23 2.316521 0.029385 0.82279
# Used Ranked model to adjust for outlier
boxplot(rank(Cl...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)
boxplot(rank(Cl...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)
# Nonparametric Wilcoxon test
wilcox.test(rank(Cl...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(Cl...mmol.L.) by Pregnant.Or.Atresia
## W = 80, p-value = 0.6463
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(rank(Cl...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(Cl...mmol.L.) by Ambient.Or.OAH
## W = 53.5, p-value = 0.2973
## alternative hypothesis: true location shift is not equal to 0
Parametric Results: Ambient Fecundity No Atresia Sample
Interaction not significant (Pr(>F) = 0.678), additive model used.
Additive summary results:
No significant differences between parturition and treatment groups.
Parametric post Hoc: Tukey HSD with additive model
No groups appear significantly different from eachother…
Evidence:
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Met, W = 0.97188, p-value = 0.6318
Homoscedasticity: Not met (for treatment group only)
Evidence:
Parturition group: Bartlett’s K-squared = 0.76057, df = 1, p-value = 0.3832
Treatment group: Bartlett’s K-squared = 3.6573, df = 1, p-value = 0.05582
Non-Parametric Results:
Kruskal-Wallis rank sum test: - No significant difference in chloride detected between factor groups. - Evidence:
Nonparametric Post-Hoc: Dunn’s test - No significant difference detected across factor groups. - Evidence: - Parturition group: pval = 0.6291 - Treatment group: pval = 0.2850
Similarities:
# chloride
# Primary Samples
# Regression Model
# p-value significant if (a < 0.1 or a < 0.05)
chloride_primary_lm_int <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(chloride_primary_lm_int)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.462 -2.050 -0.200 1.654 7.538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 158.85288 0.89826 176.845 <2e-16
## Pregnant.Or.Atresia.L -0.03128 1.27033 -0.025 0.981
## Ambient.Or.OAH.L 1.05658 1.27033 0.832 0.414
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.75577 1.79652 -0.421 0.678
##
## (Intercept) ***
## Pregnant.Or.Atresia.L
## Ambient.Or.OAH.L
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.783 on 24 degrees of freedom
## Multiple R-squared: 0.03231, Adjusted R-squared: -0.08865
## F-statistic: 0.2671 on 3 and 24 DF, p-value: 0.8484
chloride_primary_lm_add <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(chloride_primary_lm_add)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3326 -2.0528 -0.4339 1.8667 7.6674
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 158.7473 0.8482 187.160 <2e-16 ***
## Pregnant.Or.Atresia.L 0.2638 1.0415 0.253 0.802
## Ambient.Or.OAH.L 0.8503 1.1525 0.738 0.467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.721 on 25 degrees of freedom
## Multiple R-squared: 0.02517, Adjusted R-squared: -0.05281
## F-statistic: 0.3228 on 2 and 25 DF, p-value: 0.7271
# chloride_primary_lm_res_plot <- simulateResiduals(fittedModel = chloride_primary_lm_int)
# plot(chloride_primary_lm_res_plot)
chloride_primary_lm_res_plot <- simulateResiduals(fittedModel = chloride_primary_lm_add)
plot(chloride_primary_lm_res_plot) # quantile deviations detected
visualize(chloride_primary_lm_add, plot = "model")
# Check Dispersion or Outliers
plotResiduals(chloride_primary_lm_add)
testDispersion(chloride_primary_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# Transform data: square root
chloride_primary_sqrt <- Primary_Samples %>%
mutate(Cl...mmol.L. = sqrt(Cl...mmol.L.))
# Rerun lm model with square root transformed data
chloride_primary_sqrt_lm_int <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = chloride_primary_sqrt)
summary(chloride_primary_sqrt_lm_int) # still not significant
##
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = chloride_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.298717 -0.081125 -0.007769 0.066405 0.297176
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.603101 0.035643 353.591 <2e-16
## Pregnant.Or.Atresia.L -0.001406 0.050407 -0.028 0.978
## Ambient.Or.OAH.L 0.042392 0.050407 0.841 0.409
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.029585 0.071286 -0.415 0.682
##
## (Intercept) ***
## Pregnant.Or.Atresia.L
## Ambient.Or.OAH.L
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1501 on 24 degrees of freedom
## Multiple R-squared: 0.03267, Adjusted R-squared: -0.08825
## F-statistic: 0.2702 on 3 and 24 DF, p-value: 0.8462
chloride_primary_sqrt_lm_add <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = chloride_primary_sqrt)
summary(chloride_primary_sqrt_lm_add)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = chloride_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29367 -0.08112 -0.01648 0.07480 0.30222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.59897 0.03365 374.377 <2e-16 ***
## Pregnant.Or.Atresia.L 0.01015 0.04132 0.246 0.808
## Ambient.Or.OAH.L 0.03432 0.04573 0.751 0.460
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1476 on 25 degrees of freedom
## Multiple R-squared: 0.02573, Adjusted R-squared: -0.05221
## F-statistic: 0.3301 on 2 and 25 DF, p-value: 0.7219
# Check hows it
chloride_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = chloride_primary_sqrt_lm_add)
plot(chloride_primary_sqrt_lm_res_plot) # KS still significant
plotResiduals(chloride_primary_sqrt_lm_add)
testDispersion(chloride_primary_sqrt_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
chloride_primary_lm_int_res <- residuals(chloride_primary_lm_int)
chloride_primary_lm_add_res <- residuals(chloride_primary_lm_add)
# Normality test
shapiro.test(chloride_primary_lm_int_res) # normal
##
## Shapiro-Wilk normality test
##
## data: chloride_primary_lm_int_res
## W = 0.96898, p-value = 0.5534
shapiro.test(chloride_primary_lm_add_res) # normal
##
## Shapiro-Wilk normality test
##
## data: chloride_primary_lm_add_res
## W = 0.96953, p-value = 0.5681
# Scedasticity test:
# bartlett's test parametric
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.76057, df = 1, p-value = 0.3832
bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6573, df = 1, p-value = 0.05582
# Levene's test if nonparametric
leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.2426 0.3162
## 24
leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1495 0.7021
## 26
leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.9761 0.09637 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Non-parametric regression test???:
Results: Primary Sample (With both treatments and parturition groups)
Linear Model: Additive lm param ~ parturition + treatment - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.02517, Adjusted R-squared: -0.05281 F-statistic: 0.3228 on 2 and 25 DF, p-value: 0.7271
Warning: DHARMA throws error about quantile deviations (dispersion?).
Correlation: Using sqrt transformed chloride values (additive model) - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: Multiple R-squared: 0.02573, Adjusted R-squared: -0.05221 F-statistic: 0.3301 on 2 and 25 DF, p-value: 0.7219
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality (of residuals) : met
Homoscedasticity: Not Met (for treatment group)
Evidence:
Parturition group: Bartlett’s K-squared = 0.76057, df = 1, p-value = 0.3832
Treatment group: Bartlett’s K-squared = 3.6573, df = 1, p-value = 0.05582
Non-parametric Tests:
nonparametric normality test (Levene’s test): Met, sorta for treatment group
Nonparametric regression stat test???
# chloride
# Primary Samples
# Regression Model pt 2
# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
chloride_primary_lm1 <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
summary(chloride_primary_lm1)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.307152 -0.087238 0.003367 0.068150 0.288741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.58630 0.02887 436.038 <2e-16 ***
## Pregnant.Or.Atresia.L 0.01282 0.04082 0.314 0.756
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1464 on 26 degrees of freedom
## Multiple R-squared: 0.003776, Adjusted R-squared: -0.03454
## F-statistic: 0.09856 on 1 and 26 DF, p-value: 0.7561
chloride_primary_lm2 <- lm(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
summary(chloride_primary_lm2)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28820 -0.08641 -0.01170 0.07993 0.30769
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.60136 0.03163 398.457 <2e-16 ***
## Ambient.Or.OAH.L 0.03528 0.04473 0.789 0.437
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1449 on 26 degrees of freedom
## Multiple R-squared: 0.02338, Adjusted R-squared: -0.01418
## F-statistic: 0.6224 on 1 and 26 DF, p-value: 0.4373
chloride_primary_lm_res_plot1 <- simulateResiduals(fittedModel = chloride_primary_lm1)
plot(chloride_primary_lm_res_plot1)
chloride_primary_lm_res_plot2 <- simulateResiduals(fittedModel = chloride_primary_lm2)
plot(chloride_primary_lm_res_plot2)
visualize(chloride_primary_lm1, plot = "model")
visualize(chloride_primary_lm2, plot = "model")
# For shapiro test
chloride_primary_lm_res1 <- residuals(chloride_primary_lm1)
chloride_primary_lm_res2 <- residuals(chloride_primary_lm2)
# Normality test
shapiro.test(chloride_primary_lm_res1) # Normal
##
## Shapiro-Wilk normality test
##
## data: chloride_primary_lm_res1
## W = 0.97367, p-value = 0.6816
shapiro.test(chloride_primary_lm_res2) # Normal
##
## Shapiro-Wilk normality test
##
## data: chloride_primary_lm_res2
## W = 0.96785, p-value = 0.5243
# Scedasticity test: On sqrt transformed data
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.75694, df = 1, p-value = 0.3843
bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6841, df = 1, p-value = 0.05494
# Nonparametric variance test:
leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1473 0.7042
## 26
leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.0146 0.09436 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Results: Primary Sample Regression Analysis pt2
Linear Model: Individual test correlation between param and treatment, and param and parturition outcome
Correlation: Using sqrt transformed sodium values - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: - Meaning:
Warning: DHARMA says all good for square root transformed lm model
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Normal
Homoscedasticity: Not Met, for sqrt transformed data, treatment group
Evidence: On sqrt transformed data
Parturition group: Bartlett’s K-squared = 0.75694, df = 1, p-value = 0.3843
Treatment group: Bartlett’s K-squared = 3.6841, df = 1, p-value = 0.05494
Nonparametric variance test (levene’s test):
Not met for treatment group: Pr(>F) = 0.09436 .
Residuals Plot (if significant)
# chloride
# Primary samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# lm for chloride & parturition scatterplot:
chloride_primary_lm_scatterplot1 <- ggplot(data = chloride_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Parturition Group", y = "Sqrt Chloride", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_primary_lm_scatterplot1)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/chloride_primary_lm_scatterplot1.pdf", plot = chloride_primary_lm_scatterplot1, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/chloride_primary_lm_scatterplot1.png", plot = chloride_primary_lm_scatterplot1, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# lm for glucose & treatment scatterplot:
chloride_primary_lm_scatterplot2 <- ggplot(data = chloride_primary_sqrt, aes(x = Ambient.Or.OAH, y = Cl...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Treatment Group", y = "Sqrt Chloride", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_primary_lm_scatterplot2)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/chloride_primary_lm_scatterplot2.pdf", plot = chloride_primary_lm_scatterplot2, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/chloride_primary_lm_scatterplot2.png", plot = chloride_primary_lm_scatterplot2, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# chloride
# lm Boxplot
chloride_primary_lm_boxplot1 <- ggplot(data = chloride_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Parturition Group", y = "Sqrt Chloride", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_primary_lm_boxplot1)
## `geom_smooth()` using formula = 'y ~ x'
chloride_primary_lm_boxplot2 <- ggplot(data = chloride_primary_sqrt, aes(x = Ambient.Or.OAH, y = Cl...mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Parturition Group", y = "Sqrt Chloride", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_primary_lm_boxplot2)
## `geom_smooth()` using formula = 'y ~ x'
Potassium Summary Stats
# K+
# Summary stats
# Primary data
Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarize(count = n(),
median = round(median(K...mmol.L.), 3),
mean = round(mean(K...mmol.L.), 3),
sd = round(sd(K...mmol.L.), 3),
cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.
## # A tibble: 4 × 7
## Pregnant.Or.Atresia Ambient.Or.OAH count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition Ambient 8 2.7 2.71 0.146 0.054
## 2 Post Parturition OAH 2 3 3 0.283 0.094
## 3 Atresia Ambient 13 2.7 2.86 0.304 0.106
## 4 Atresia OAH 5 2.7 3 0.714 0.238
Potassium Boxplot
# K+ boxplot
# Primary samples
potassium_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = K...mmol.L., fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = K...mmol.L., fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
labs(title = "Potassium", x = "Parturition Type", y = "Potassium (mmol/L)") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_primary_boxplot)
ggsave(filename = "data-images/potassium_primary_boxplot.pdf", plot = potassium_primary_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_primary_boxplot.png", plot = potassium_primary_boxplot, device = "png")
## Saving 7 x 5 in image
Data Distribution
# K+
# Parametric Assumptions
# Data Distribution
hist(Primary_Samples$K...mmol.L.)
plotNormalHistogram(Primary_Samples$K...mmol.L.)
plotNormalDensity(Primary_Samples$K...mmol.L.)
Differences:
# K+
# Primary Samples
# Interactive aov model
potassium_primary_aov_int <- aov(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_aov_int) # interaction not significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 0.109 0.10864 0.772 0.388
## Ambient.Or.OAH 1 0.177 0.17685 1.256 0.274
## Pregnant.Or.Atresia:Ambient.Or.OAH 1 0.025 0.02463 0.175 0.680
## Residuals 24 3.380 0.14081
# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
potassium_primary_aov_add <- aov(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_aov_add)
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 0.109 0.1086 0.798 0.380
## Ambient.Or.OAH 1 0.177 0.1769 1.299 0.265
## Residuals 25 3.404 0.1362
# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(potassium_primary_aov_add)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
##
## $Pregnant.Or.Atresia
## diff lwr upr p adj
## Atresia-Post Parturition 0.13 -0.1697415 0.4297415 0.3802489
##
## $Ambient.Or.OAH
## diff lwr upr p adj
## OAH-Ambient 0.1828571 -0.1488268 0.5145411 0.266967
# Normality plot
# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
potassium_primary_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Potassium Interactive Residual QQplot",
subtitle = "residuals(aov(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
x = "Potassium Theoretical Quantiles (Predicted values)",
y = "Potassium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5))
print(potassium_primary_res_qqplot)
# Normality Test
shapiro.test(Primary_Samples$K...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Primary_Samples$K...mmol.L.
## W = 0.76701, p-value = 2.949e-05
shapiro.test(residuals(lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.86428, p-value = 0.00184
# Parametric variance test:
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 5.4475, df = 1, p-value = 0.0196
bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 7.2218, df = 1, p-value = 0.007202
# Nonparametric variance test:
LeveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3073 0.2633
## 26
leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.3217 0.07989 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Nonparametric Stat test & Post-Hoc:
kruskal.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: K...mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 0.15497, df = 1, p-value = 0.6938
kruskal.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: K...mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 0.24017, df = 1, p-value = 0.6241
DunnTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Atresia-Post Parturition 1.244444 0.6938
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DunnTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## OAH-Ambient 1.714286 0.6241
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
visualize(potassium_primary_aov_add)
# Test outliers
outlierTest(potassium_primary_aov_add)
## rstudent unadjusted p-value Bonferroni p
## 17 4.648451 0.00010152 0.0028426
# Used Ranked model to adjust for outlier
boxplot(rank(K...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)
boxplot(rank(K...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)
# Nonparametric Wilcoxon test
wilcox.test(rank(K...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(K...mmol.L.) by Pregnant.Or.Atresia
## W = 82, p-value = 0.7121
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(rank(K...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(K...mmol.L.) by Ambient.Or.OAH
## W = 64.5, p-value = 0.6435
## alternative hypothesis: true location shift is not equal to 0
Parametric Results: Ambient Fecundity No Atresia Sample
Interaction not significant (Pr(>F) = 0.680), additive model used.
Additive summary results:
No significant differences between parturition and treatment groups.
Parametric post Hoc: Tukey HSD with additive model
No groups appear significantly different from eachother…
Evidence:
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality: Not Met
Homoscedasticity: Not met (for both factors)
Evidence:
Parturition group:
Treatment group:
Non-Parametric Results:
Kruskal-Wallis rank sum test: - No significant difference in chloride detected between parturition and treatment factors - Evidence:
Nonparametric Post-Hoc: Dunn’s test - No significant difference detected across factor groups. - Evidence: - Parturition group: - Treatment group:
Outlier test: Appears significant
Similarities:
# potassium
# Primary Samples
# Regression Model
# p-value significant if (a < 0.1 or a < 0.05)
potassium_primary_lm_int <- lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(potassium_primary_lm_int)
##
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50000 -0.17115 -0.03702 0.09063 1.20000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.89351 0.08909 32.477 <2e-16 ***
## Pregnant.Or.Atresia.L 0.05269 0.12600 0.418 0.680
## Ambient.Or.OAH.L 0.15060 0.12600 1.195 0.244
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.07452 0.17819 -0.418 0.680
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3753 on 24 degrees of freedom
## Multiple R-squared: 0.08405, Adjusted R-squared: -0.03044
## F-statistic: 0.7341 on 3 and 24 DF, p-value: 0.5419
potassium_primary_lm_add <- lm(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_lm_add)
##
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53305 -0.14883 -0.04099 0.06687 1.16695
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.88310 0.08412 34.272 <2e-16 ***
## Pregnant.Or.Atresia.L 0.08179 0.10329 0.792 0.436
## Ambient.Or.OAH.L 0.13026 0.11430 1.140 0.265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.369 on 25 degrees of freedom
## Multiple R-squared: 0.07738, Adjusted R-squared: 0.003568
## F-statistic: 1.048 on 2 and 25 DF, p-value: 0.3654
# potassium_primary_lm_res_plot <- simulateResiduals(fittedModel = potassium_primary_lm_int)
# plot(potassium_primary_lm_res_plot)
potassium_primary_lm_res_plot <- simulateResiduals(fittedModel = potassium_primary_lm_add)
plot(potassium_primary_lm_res_plot) # quantile deviations detected
visualize(potassium_primary_lm_add, plot = "model")
# Check Dispersion or Outliers
plotResiduals(potassium_primary_lm_add)
testDispersion(potassium_primary_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# For shapiro.test
potassium_primary_lm_add_res <- residuals(potassium_primary_lm_add)
# Normality test
shapiro.test(potassium_primary_lm_add_res) # nonparametric
##
## Shapiro-Wilk normality test
##
## data: potassium_primary_lm_add_res
## W = 0.87247, p-value = 0.002751
# Scedasticity test:
# bartlett's test parametric
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 5.4475, df = 1, p-value = 0.0196
bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 7.2218, df = 1, p-value = 0.007202
# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3073 0.2633
## 26
leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.3217 0.07989 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Results: Primary Sample (With both treatments and parturition groups)
Linear Model: Additive lm param ~ parturition + treatment - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.07738, Adjusted R-squared: 0.003568 F-statistic: 1.048 on 2 and 25 DF, p-value: 0.3654
Warning: No DHARMA errors (except a star).
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality (of residuals) : Not met
Homoscedasticity: Not Met
Evidence:
Parturition group:
Treatment group:
Non-parametric results:
nonparametric variance test (Levene’s test): - Homoscedasticity: Not met for treatment group
transform data (if necessary)
# potassium
# Primary Samples
# Transform data: square root
potassium_primary_sqrt <- Primary_Samples %>%
mutate(K...mmol.L. = sqrt(K...mmol.L.))
# Rerun lm model with square root transformed data
potassium_primary_sqrt_lm_int <- lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = potassium_primary_sqrt)
summary(potassium_primary_sqrt_lm_int) # still not significant
##
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = potassium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14196 -0.04920 -0.00974 0.02954 0.32629
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.69754 0.02503 67.819 <2e-16 ***
## Pregnant.Or.Atresia.L 0.01240 0.03540 0.350 0.729
## Ambient.Or.OAH.L 0.04180 0.03540 1.181 0.249
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.02552 0.05006 -0.510 0.615
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1054 on 24 degrees of freedom
## Multiple R-squared: 0.08134, Adjusted R-squared: -0.03349
## F-statistic: 0.7084 on 3 and 24 DF, p-value: 0.5564
potassium_primary_sqrt_lm_add <- lm(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = potassium_primary_sqrt)
summary(potassium_primary_sqrt_lm_add)
##
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = potassium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15328 -0.04199 -0.01110 0.02141 0.31497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.69398 0.02368 71.549 <2e-16 ***
## Pregnant.Or.Atresia.L 0.02237 0.02907 0.769 0.449
## Ambient.Or.OAH.L 0.03483 0.03217 1.083 0.289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1039 on 25 degrees of freedom
## Multiple R-squared: 0.0714, Adjusted R-squared: -0.002893
## F-statistic: 0.9611 on 2 and 25 DF, p-value: 0.3962
# Check hows it
potassium_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = potassium_primary_sqrt_lm_add)
plot(potassium_primary_sqrt_lm_res_plot) # KS still significant
plotResiduals(potassium_primary_sqrt_lm_add)
testDispersion(potassium_primary_sqrt_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
potassium_primary_lm_add_sqrt_res <- residuals(potassium_primary_sqrt_lm_add)
# Normality test
shapiro.test(potassium_primary_lm_add_sqrt_res) # nonparametric
##
## Shapiro-Wilk normality test
##
## data: potassium_primary_lm_add_sqrt_res
## W = 0.89113, p-value = 0.007134
# Scedasticity test:
# bartlett's test parametric
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = potassium_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 4.737, df = 1, p-value = 0.02952
bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = potassium_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 6.3646, df = 1, p-value = 0.01164
# Levene's test if nonparametric
leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = potassium_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.4427 0.2551
## 24
leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = potassium_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3218 0.2607
## 26
leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = potassium_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.347 0.07882 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Non-parametric regression test???:
Correlation: Using sqrt transformed chloride values (additive model) - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence:
Parametric Assumptions: Across Sample comparisons…
Sample size (n > 30): Not Met
Normality (of residuals) : Not met
Homoscedasticity: Not Met (for both groups)
Evidence:
Parturition group:
Treatment group:
Non-parametric Tests:
nonparametric normality test (Levene’s test): Not really met for treatment group
Nonparametric regression stat test???
# chloride
# Primary Samples
# Regression Model pt 2
# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
chloride_primary_lm1 <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
summary(chloride_primary_lm1)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.307152 -0.087238 0.003367 0.068150 0.288741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.58630 0.02887 436.038 <2e-16 ***
## Pregnant.Or.Atresia.L 0.01282 0.04082 0.314 0.756
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1464 on 26 degrees of freedom
## Multiple R-squared: 0.003776, Adjusted R-squared: -0.03454
## F-statistic: 0.09856 on 1 and 26 DF, p-value: 0.7561
chloride_primary_lm2 <- lm(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
summary(chloride_primary_lm2)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28820 -0.08641 -0.01170 0.07993 0.30769
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.60136 0.03163 398.457 <2e-16 ***
## Ambient.Or.OAH.L 0.03528 0.04473 0.789 0.437
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1449 on 26 degrees of freedom
## Multiple R-squared: 0.02338, Adjusted R-squared: -0.01418
## F-statistic: 0.6224 on 1 and 26 DF, p-value: 0.4373
chloride_primary_lm_res_plot1 <- simulateResiduals(fittedModel = chloride_primary_lm1)
plot(chloride_primary_lm_res_plot1)
chloride_primary_lm_res_plot2 <- simulateResiduals(fittedModel = chloride_primary_lm2)
plot(chloride_primary_lm_res_plot2)
visualize(chloride_primary_lm1, plot = "model")
visualize(chloride_primary_lm2, plot = "model")
# For shapiro test
chloride_primary_lm_res1 <- residuals(chloride_primary_lm1)
chloride_primary_lm_res2 <- residuals(chloride_primary_lm2)
# Normality test
shapiro.test(chloride_primary_lm_res1) # Normal
##
## Shapiro-Wilk normality test
##
## data: chloride_primary_lm_res1
## W = 0.97367, p-value = 0.6816
shapiro.test(chloride_primary_lm_res2) # Normal
##
## Shapiro-Wilk normality test
##
## data: chloride_primary_lm_res2
## W = 0.96785, p-value = 0.5243
# Scedasticity test: On sqrt transformed data
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.75694, df = 1, p-value = 0.3843
bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6841, df = 1, p-value = 0.05494
# Nonparametric variance test:
leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1473 0.7042
## 26
leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.0146 0.09436 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# potassium
# Primary Samples
# Regression Model pt 2: Not using sqrt model
# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
potassium_primary_lm1 <- lm(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
summary(potassium_primary_lm1)
##
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4000 -0.2000 -0.0700 0.0725 1.3000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.83500 0.07319 38.737 <2e-16 ***
## Pregnant.Or.Atresia.L 0.09192 0.10350 0.888 0.383
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3711 on 26 degrees of freedom
## Multiple R-squared: 0.02945, Adjusted R-squared: -0.007884
## F-statistic: 0.7888 on 1 and 26 DF, p-value: 0.3826
potassium_primary_lm2 <- lm(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_lm2)
##
## Call:
## lm(formula = K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50000 -0.20119 -0.10476 0.09643 1.20000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.90238 0.07994 36.305 <2e-16 ***
## Ambient.Or.OAH.L 0.13805 0.11306 1.221 0.233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3664 on 26 degrees of freedom
## Multiple R-squared: 0.05424, Adjusted R-squared: 0.01786
## F-statistic: 1.491 on 1 and 26 DF, p-value: 0.233
potassium_primary_lm_res_plot1 <- simulateResiduals(fittedModel = potassium_primary_lm1)
plot(potassium_primary_lm_res_plot1)
potassium_primary_lm_res_plot2 <- simulateResiduals(fittedModel = potassium_primary_lm2)
plot(potassium_primary_lm_res_plot2)
visualize(potassium_primary_lm1, plot = "model")
visualize(potassium_primary_lm2, plot = "model")
# For shapiro test
potassium_primary_lm_res1 <- residuals(potassium_primary_lm1)
potassium_primary_lm_res2 <- residuals(potassium_primary_lm2)
# Normality test
shapiro.test(potassium_primary_lm_res1) # Normal
##
## Shapiro-Wilk normality test
##
## data: potassium_primary_lm_res1
## W = 0.81788, p-value = 0.0002244
shapiro.test(potassium_primary_lm_res2) # Normal
##
## Shapiro-Wilk normality test
##
## data: potassium_primary_lm_res2
## W = 0.84231, p-value = 0.0006565
# Scedasticity test: On sqrt transformed data
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 5.4475, df = 1, p-value = 0.0196
bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 7.2218, df = 1, p-value = 0.007202
# Nonparametric variance test:
leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3073 0.2633
## 26
leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.3217 0.07989 .
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residuals Plot (if significant)
# potassium
# Primary samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# lm for chloride & parturition scatterplot:
potassium_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = K...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Parturition Group", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_primary_lm_scatterplot1)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/potassium_primary_lm_scatterplot1.pdf", plot = potassium_primary_lm_scatterplot1, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/potassium_primary_lm_scatterplot1.png", plot = potassium_primary_lm_scatterplot1, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# lm for glucose & treatment scatterplot:
potassium_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = K...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Treatment Group", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_primary_lm_scatterplot2)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/potassium_primary_lm_scatterplot2.pdf", plot = potassium_primary_lm_scatterplot2, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/potassium_primary_lm_scatterplot2.png", plot = potassium_primary_lm_scatterplot2, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# K+
# lm Boxplot
potassium_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = K...mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Parturition Group", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_primary_lm_boxplot1)
## `geom_smooth()` using formula = 'y ~ x'
potassium_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = K...mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Parturition Group", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_primary_lm_boxplot2)
## `geom_smooth()` using formula = 'y ~ x'
Calcium Summary Stats
# Ca+2
# Summary stats
# Primary data
Primary_Samples %>%
group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
summarize(count = n(),
median = round(median(Ca....mmol.L.), 3),
mean = round(mean(Ca....mmol.L.), 3),
sd = round(sd(Ca....mmol.L.), 3),
cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.
## # A tibble: 4 × 7
## Pregnant.Or.Atresia Ambient.Or.OAH count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition Ambient 8 1.43 1.43 0.041 0.029
## 2 Post Parturition OAH 2 1.25 1.25 0.071 0.057
## 3 Atresia Ambient 13 1.24 1.27 0.114 0.09
## 4 Atresia OAH 5 1.16 1.16 0.031 0.027
Calcium Boxplot
# Ca+2 boxplot
# Primary samples
calcium_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L., fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L., fill = Ambient.Or.OAH)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
labs(title = "Calcium", x = "Parturition Type", y = "Calcium (mmol/L)") +
guides(fill = guide_legend((title = "Treatment Type"))) +
scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_primary_boxplot)
ggsave(filename = "data-images/calcium_primary_boxplot.pdf", plot = calcium_primary_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_primary_boxplot.png", plot = calcium_primary_boxplot, device = "png")
## Saving 7 x 5 in image
Data Distribution
# Ca+2
# Parametric Assumptions
# Data Distribution
hist(Primary_Samples$Ca....mmol.L.)
plotNormalHistogram(Primary_Samples$Ca....mmol.L.)
plotNormalDensity(Primary_Samples$Ca....mmol.L.)
Differences:
# Ca+2
# Primary Samples
# Interactive aov model
calcium_primary_aov_int <- aov(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_aov_int) # interaction not significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 0.15667 0.15667 21.391 0.000108 ***
## Ambient.Or.OAH 1 0.08990 0.08990 12.274 0.001827 **
## Pregnant.Or.Atresia:Ambient.Or.OAH 1 0.00575 0.00575 0.785 0.384354
## Residuals 24 0.17578 0.00732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
calcium_primary_aov_add <- aov(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_aov_add)
## Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia 1 0.1567 0.15667 21.58 9.35e-05 ***
## Ambient.Or.OAH 1 0.0899 0.08990 12.38 0.00169 **
## Residuals 25 0.1815 0.00726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(calcium_primary_aov_add)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
##
## $Pregnant.Or.Atresia
## diff lwr upr p adj
## Atresia-Post Parturition -0.1561111 -0.2253289 -0.08689334 9.35e-05
##
## $Ambient.Or.OAH
## diff lwr upr p adj
## OAH-Ambient -0.1303704 -0.2069644 -0.05377629 0.0017412
# Normality plot
# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
calcium_primary_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Calcium Interactive Residual QQplot",
subtitle = "residuals(aov(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
x = "Calcium Theoretical Quantiles (Predicted values)",
y = "Calcium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5))
print(calcium_primary_res_qqplot)
# Normality Test
shapiro.test(Primary_Samples$Ca....mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Primary_Samples$Ca....mmol.L.
## W = 0.93132, p-value = 0.06655
shapiro.test(residuals(lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.91748, p-value = 0.03007
# Parametric variance test:
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.49921, df = 1, p-value = 0.4798
bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.5462, df = 1, p-value = 0.05968
# Nonparametric variance test:
LeveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2268 0.6379
## 26
leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 6.801 0.0149 *
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Nonparametric Stat test & Post-Hoc:
kruskal.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Ca....mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 9.7286, df = 1, p-value = 0.001814
kruskal.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Ca....mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 8.0703, df = 1, p-value = 0.0045
DunnTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Atresia-Post Parturition -10.11111 0.0018 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
DunnTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## OAH-Ambient -10.19048 0.0045 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)
boxplot(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)
# Test outliers
outlierTest(calcium_primary_aov_add)
## rstudent unadjusted p-value Bonferroni p
## 16 3.869496 0.00073227 0.020504
# Used Ranked model to adjust for outlier
boxplot(rank(Ca....mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)
boxplot(rank(Ca....mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)
# Nonparametric Wilcoxon test
wilcox.test(rank(Ca....mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, alternative = c("two.sided"))
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(Ca....mmol.L.) by Pregnant.Or.Atresia
## W = 155, p-value = 0.001968
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(rank(Ca....mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, alternative = c("two.sided"))
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: rank(Ca....mmol.L.) by Ambient.Or.OAH
## W = 127, p-value = 0.004889
## alternative hypothesis: true location shift is not equal to 0
Interaction not significant, additive model used.
Results: summary aov additive model - Post parturition group significantly different from atresia group (Pr(>F) = 9.35e-05 *) - Ambient treatment significantly different from OAH treatment (Pr(>F) = 0.00169 )
Parametric Post Hoc: Tukey HSD test - Post parturition patients display a significant increase in calcium levels compared to atresia patients (p adj = 9.35e-05) - Ambient treatment patients display a significant increase in calcium levels compared to OAH treatment patients (p adj = 0.0017412)
Parametric Assumptions:
Normality: Not met - W = 0.93132, p-value = 0.06655
Parametric Scedasticity test (Bartlett’s test): - Not met for treatment groups
Nonparametric variance test: Levene’s test: - Not met for treatment group
Nonparametric Stat test results: Kruskall Wallis test - Post parturition patients display a significant increase in calcium levels compared to atresia patients (Kruskal-Wallis chi-squared = 9.7286, df = 1, p-value = 0.001814) - Ambient treatment patients display a significant increase in calcium levels compared to OAH treatment patients (Kruskal-Wallis chi-squared = 8.0703, df = 1, p-value = 0.0045)
Nonparametric Post Hoc: Dunn’s Test - nparametric Stat test results: Kruskall Wallis test - Post parturition patients display a significant increase in calcium levels compared to atresia patients (OAH-Ambient: mean.rank.diff = -10.19048, pval = 0.0018 ) - Ambient treatment patients display a significant increase in calcium levels compared to OAH treatment patients (OAH-Ambient: mean.rank.diff = -10.19048, pval = 0.0045 )
Similarities:
# calcium
# Primary Samples
# Regression Model
# p-value significant if (a < 0.1 or a < 0.05)
calcium_primary_lm_int <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(calcium_primary_lm_int)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14923 -0.05000 -0.00524 0.03219 0.26077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.27762 0.02032 62.878 < 2e-16 ***
## Pregnant.Or.Atresia.L -0.08910 0.02874 -3.101 0.00488 **
## Ambient.Or.OAH.L -0.10270 0.02874 -3.574 0.00153 **
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L 0.03601 0.04064 0.886 0.38435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08558 on 24 degrees of freedom
## Multiple R-squared: 0.5894, Adjusted R-squared: 0.5381
## F-statistic: 11.48 on 3 and 24 DF, p-value: 7.291e-05
calcium_primary_lm_add <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_lm_add)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.155373 -0.042295 -0.003321 0.038470 0.254627
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.28265 0.01943 66.026 < 2e-16 ***
## Pregnant.Or.Atresia.L -0.10316 0.02385 -4.325 0.000214 ***
## Ambient.Or.OAH.L -0.09287 0.02640 -3.519 0.001685 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08521 on 25 degrees of freedom
## Multiple R-squared: 0.576, Adjusted R-squared: 0.542
## F-statistic: 16.98 on 2 and 25 DF, p-value: 2.201e-05
# calcium_primary_lm_res_plot <- simulateResiduals(fittedModel = calcium_primary_lm_int)
# plot(calcium_primary_lm_res_plot)
calcium_primary_lm_res_plot <- simulateResiduals(fittedModel = calcium_primary_lm_add)
plot(calcium_primary_lm_res_plot) # quantile deviations detected
visualize(calcium_primary_lm_add, plot = "model")
# Check Dispersion or Outliers
plotResiduals(calcium_primary_lm_add)
testDispersion(calcium_primary_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# For shapiro.test
calcium_primary_lm_add_res <- residuals(calcium_primary_lm_add)
# Normality test
shapiro.test(calcium_primary_lm_add_res) # nonparametric
##
## Shapiro-Wilk normality test
##
## data: calcium_primary_lm_add_res
## W = 0.94439, p-value = 0.1429
# Scedasticity test:
# bartlett's test parametric
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.49921, df = 1, p-value = 0.4798
bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.5462, df = 1, p-value = 0.05968
# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2268 0.6379
## 26
leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 6.801 0.0149 *
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Interaction not significant, additive model used.
lm summary results: additive model - Significant relationship detected. - Multiple R-squared: 0.576, Adjusted R-squared: 0.542 F-statistic: 16.98 on 2 and 25 DF, p-value: 2.201e-05
Parametric Assumptions:
Normality: Not met
Parametric Scedasticity test (Bartlett’s test): - Not met for treatment groups
Nonparametric variance test: Levene’s test: - Not met for treatment group
# Ca+2
# Primary Samples
# Regression Model pt 2: Not using sqrt model
# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
calcium_primary_lm1 <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
summary(calcium_primary_lm1)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.195000 -0.063889 -0.006944 0.040000 0.291111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.31694 0.02015 65.360 < 2e-16 ***
## Pregnant.Or.Atresia.L -0.11039 0.02849 -3.874 0.000649 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1022 on 26 degrees of freedom
## Multiple R-squared: 0.366, Adjusted R-squared: 0.3416
## F-statistic: 15.01 on 1 and 26 DF, p-value: 0.000649
calcium_primary_lm2 <- lm(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_lm2)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.210952 -0.083452 0.001667 0.091548 0.199048
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.25833 0.02411 52.194 < 2e-16 ***
## Ambient.Or.OAH.L -0.10270 0.03409 -3.012 0.00571 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1105 on 26 degrees of freedom
## Multiple R-squared: 0.2587, Adjusted R-squared: 0.2302
## F-statistic: 9.073 on 1 and 26 DF, p-value: 0.005714
calcium_primary_lm_res_plot1 <- simulateResiduals(fittedModel = calcium_primary_lm1)
plot(calcium_primary_lm_res_plot1)
calcium_primary_lm_res_plot2 <- simulateResiduals(fittedModel = calcium_primary_lm2)
plot(calcium_primary_lm_res_plot2)
visualize(calcium_primary_lm1, plot = "model")
visualize(calcium_primary_lm2, plot = "model")
# Check Dispersion or Outliers
plotResiduals(calcium_primary_lm1)
testDispersion(calcium_primary_lm1)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97355, p-value = 0.928
## alternative hypothesis: two.sided
plotResiduals(calcium_primary_lm2)
testDispersion(calcium_primary_lm2)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97355, p-value = 0.928
## alternative hypothesis: two.sided
# For shapiro test
calcium_primary_lm_res1 <- residuals(potassium_primary_lm1)
calcium_primary_lm_res2 <- residuals(potassium_primary_lm2)
# Normality test
shapiro.test(calcium_primary_lm_res1) # Normal
##
## Shapiro-Wilk normality test
##
## data: calcium_primary_lm_res1
## W = 0.81788, p-value = 0.0002244
shapiro.test(calcium_primary_lm_res2) # Normal
##
## Shapiro-Wilk normality test
##
## data: calcium_primary_lm_res2
## W = 0.84231, p-value = 0.0006565
# Scedasticity test: On sqrt transformed data
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.49921, df = 1, p-value = 0.4798
bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.5462, df = 1, p-value = 0.05968
# Nonparametric variance test:
leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2268 0.6379
## 26
leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 6.801 0.0149 *
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Results: - A significant correlation detected, increasing calcium levels associated with post-parturition patients opposed to atresia patients.
Warning DHAMA error thrown for treatment group: Levene test for homogeneity of variance significant
Normality: Not met (for both factor groups)
Parametric Equal Variance test: Bartlett’s test - Not met for treatment group
Nonparametric variance test: Levene’s test - Not met for treatment group
Residuals Plot (if significant)
# Ca+2
# Primary samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# lm for chloride & parturition scatterplot:
calcium_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Parturition Group", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_primary_lm_scatterplot1)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/calcium_primary_lm_scatterplot1.pdf", plot = calcium_primary_lm_scatterplot1, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/calcium_primary_lm_scatterplot1.png", plot = calcium_primary_lm_scatterplot1, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# lm for glucose & treatment scatterplot:
calcium_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Ca....mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Treatment Group", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_primary_lm_scatterplot2)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/calcium_primary_lm_scatterplot2.pdf", plot = calcium_primary_lm_scatterplot2, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/calcium_primary_lm_scatterplot2.png", plot = calcium_primary_lm_scatterplot2, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# Ca+2
# lm Boxplot
calcium_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Parturition Group", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_primary_lm_boxplot1)
## `geom_smooth()` using formula = 'y ~ x'
calcium_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Ca....mmol.L.)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Parturition Group", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_primary_lm_boxplot2)
## `geom_smooth()` using formula = 'y ~ x'
transform data (if necessary)
# Ca+2
# Primary Samples
# Transform data: square root
calcium_primary_sqrt <- Primary_Samples %>%
mutate(Ca....mmol.L. = sqrt(Ca....mmol.L.))
# Rerun lm model with square root transformed data
calcium_primary_sqrt_lm_int <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = calcium_primary_sqrt)
summary(calcium_primary_sqrt_lm_int) # still not significant
##
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH,
## data = calcium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.067294 -0.021721 -0.001756 0.014485 0.111337
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.129150 0.008855 127.513 < 2e-16
## Pregnant.Or.Atresia.L -0.039421 0.012523 -3.148 0.00436
## Ambient.Or.OAH.L -0.044925 0.012523 -3.587 0.00148
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L 0.014896 0.017710 0.841 0.40860
##
## (Intercept) ***
## Pregnant.Or.Atresia.L **
## Ambient.Or.OAH.L **
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0373 on 24 degrees of freedom
## Multiple R-squared: 0.5924, Adjusted R-squared: 0.5415
## F-statistic: 11.63 on 3 and 24 DF, p-value: 6.689e-05
calcium_primary_sqrt_lm_add <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = calcium_primary_sqrt)
summary(calcium_primary_sqrt_lm_add)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH,
## data = calcium_primary_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.069835 -0.017804 -0.000962 0.017027 0.108797
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.131231 0.008453 133.829 < 2e-16 ***
## Pregnant.Or.Atresia.L -0.045238 0.010379 -4.359 0.000196 ***
## Ambient.Or.OAH.L -0.040860 0.011485 -3.558 0.001527 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03708 on 25 degrees of freedom
## Multiple R-squared: 0.5804, Adjusted R-squared: 0.5468
## F-statistic: 17.29 on 2 and 25 DF, p-value: 1.93e-05
# Check hows it
calcium_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = calcium_primary_sqrt_lm_add)
plot(calcium_primary_sqrt_lm_res_plot) # KS still significant
plotResiduals(calcium_primary_sqrt_lm_add)
testDispersion(calcium_primary_sqrt_lm_add)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
calcium_primary_lm_add_sqrt_res <- residuals(calcium_primary_sqrt_lm_add)
# Normality test
shapiro.test(calcium_primary_lm_add_sqrt_res) # nonparametric
##
## Shapiro-Wilk normality test
##
## data: calcium_primary_lm_add_sqrt_res
## W = 0.94977, p-value = 0.1955
# Scedasticity test:
# bartlett's test parametric
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = calcium_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.57932, df = 1, p-value = 0.4466
bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = calcium_primary_sqrt)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.1949, df = 1, p-value = 0.07387
# Levene's test if nonparametric
leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = calcium_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.7959 0.1749
## 24
leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = calcium_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.327 0.5723
## 26
leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = calcium_primary_sqrt)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 6.1011 0.02039 *
## 26
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Non-parametric regression test???: